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We review the main properties of a supersolid. We describe first the macroscopic equation that 
satisfies a supersolid based on general arguments and symmetries and show that such sohds might 
exhibit simultaneously or independently both elastic behavior and superfluidity. We then explain 
why a supersolid state should exist for solids at very low temperature but with a very small superfiuid 
fraction. Finally, we propose a mean-field model, based on the Gross-Pitaevskii equation, which 
presents the general properties expected for a supersolid and should therefore provide a consistent 
framework to study its dynamical properties. 
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I. INTRODUCTION 

The striking observation by Kim and Chan in 2004 of a rapid drop in the rotational inertial of solid Helium below 
O.lKelvin [Tl[2] has revived the question whether the solids could exhibit a new state of the matter at low temperature, 
the so-called "supersolidity". In these experiments, the supersolidity appears as a decoupling between a "superfluid" 
and a "normal" part of the solid: everything happens as if a small fraction of the total mass of solid Helium (of 
the order of a few percent) is at rest while a rotation is imposed to the sample! In fact, since the liquefaction of 
Helium at T = 4.22 Kelvin in 1908 by Kammerlingh Onnes, the physics of Helium has accompanied many important 
scientific discoveries in the last century. Among them, the remarkable property of superfluidity was observed in 1938 
for Helium 4 (the bosonic and most abundant isotope of Helium, noted also ''He) by P.L. Kapitza OH]. Below the 
A-point separating the two Helium liquid phases noted Helium I (above 2.2 Kelvin) and Helium II (below) [5l[6], the 
resistance to flow of Helium through thin capillaries drops suddenly so that it can be considered a zero- viscosity liquid. 
Kapitza named this effect superfluidity by analogy with the superconductivity although Helium 4 is a boson while 
superconductivity concerns fermions, the electrons. A slightly later Allen and Jones [7j exhibited the fountain effect, 
a macroscopic striking consequence of the superfluidity. It is generally believed that superfluidity is related to the 
Bose- Einstein (B-E) condensation [5] [3] and both the superconductivity and the superfluidity of Helium 3 (a fermion) 
are explained by the formation of Cooper pairs which are bosons. However, no direct connection between superfluidity 
and Bose-Einstein condensation is invoked in the classical theory of L.D. Landau [101 [H] where the dynamics of a 
so-called "coherent state" is deduced within the framework of general quantum many body theory, with no reference 
to the quantum statistics, fermionic or bosonic. Let us remark that in this respect the Bose-Einstein condensate is 
not superfluid, while it becomes superfluid when interactions are taken into account [T^. Landau is led to introduce 
then the concept of a two-fluids model: the liquid Helium below the A-point can be decomposed between a superfluid 
part which has zero viscosity and a normal part which consists of the thermal excitations of the liquid. 

On the other hand, the solidification of Helium being achieved in 1926 by Keesom, the ability for solid Helium to 
exhibit B-E condensation has been naturally questioned !l3'. Later on, in a seminal paper where they link Bose- 
Einstein condensation (and thus in their spirit superfluidity) to off diagonal long range order (ODLRO), Penrose 
and Onsager [T3] have investigated the possibility of a superfluid-like behavior in a solid at low temperature, finally 
for concluding to the inconsistency for a perfect crystal of such a "supersolid" state: the short range solid order 
is incompatible with a long range order. Further works by Andreev and Lifshitz [TS], Reatto [H] and Chester [IT] 
have suggested that supersolids, if they exist, might be a kind of Bose-Einstein condensation of defects, vacancies, 
interstitials, etc. that display a coherent state and can realize a matter flow through the crystal. Following the theory 
by Andreev and Lifshitz I15j such a superflow was expected under forces due to an imposed pressure gradient, likewise 
what had been seen by Kapitza in liquid Helium. The very first experiment by Andreev and collaborators in [18] did 
not show superflow generated in this way: they put a small sphere inside a bucket filled of solid helium, expecting it 
to fall down across the solid under its own weight in the field of gravity. Similar experiments have shown as well that 
solid Helium does not flow under an imposed pressure difference [19, 20J. 

In likely the most important theoretical paper on supersolids since Penrose and Onsager's first publication, Leggett 
suggested l2l! that non classical rotational inertia (NCRI) may be a test for supersolid. NCRI consists in an observable 
drop of the effective moment of inertia of a substance in a container of given geometry, like annular, cylindrical, cubic 
or porous vycor. This moment of inertia is measured sensitively by the frequency of resonance of a torsion pendulum, 
where the mass is a piece of supersolid. Such a drop is caused by the fact that the superfluid component of the 
supersolid does not follow the rotational motion of the rest of the solid, although the details depend on the geometry 
of the solid. Such effect is of course inspired by the classical measurement of superfluid density of liquid Helium by 
Androniksahvili [22^ 23, . In a NCRI experiment the effective moment of inertia (which is always less than the one of a 
rigid body, the difference being due to the nonconvected superfluid fraction) is determined for different temperatures 
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thanks to a torsional resonator. Once geometrical effects are taken into account, measurements of NCRI give access 
to the superfluid density which should be a universal function of the thermodynamical parameters, pressure and 
temperature. It happens that, experimentally, and contrary to what happens with superfluid liquid helium, the 
superfluid fraction of a supersolid is not a function of the thermodynamic parameters only, a point we discuss at 
length below. 

Although non classical rotational inertia was expected to be a fundamental characteristic of supersolids, experiments 
by Reppy and collaborators 24 in the early 80's failed to exhibit any measurable NCRI in solid helium at cryogenic 
temperatures. In fact, over the last 40 years a more or less ongoing research activity has failed to put in evidence 
in a convincing way the supersolid state of matter |2S], either by NCRI or through eventual changes in the modes 
of propagation of disturbances like sound waves brought by the superfluid component |26j . But recently Chan and 
Kim [1] [21 [371 [2H] have claimed the observation of NCRI in solid Helium 4 below 0.1 — 0.2K, something which gave 
new life to the subject. Measuring the resonant frequency of a torsional oscillator filled of solid Helium, they observe 
a rapid drop of the rotational moment of inertia of the order of a few percent when decreasing the temperature below 
O.IK. The successful difference between Chan and Kim's experiments and the previous attempts lied in particular in 
the very small amplitude of the oscillations achieved to measure the momentum of inertia. To date, at least five other 
groups have confirmed Kim and Chan findings |29H331 [55] but showing important variations in the superfluid fraction 
with the experimental configuration as already noticed by Kim and Chan. In particular, crystal annealing is shown 
to lower dramatically the superfluid fraction with a strong dependence on the ratio between the surface and the bulk 
of the sample (something quantified by the disorder in the solid crystal) up to a variation of almost three orders of 
magnitude for the NCRI fraction (NCRIF) [29, 31]. In addition, no evidence of superflow has been noticed in solid 
He"* submitted to a localized pressure jump [35n37] in the conditions where NCRI exists. For these reasons, the status 
of this NCRI effect as a proof of a new "supersolid" state of matter has been often questioned and remains a subject 
of scientific debates. However, we would like to emphasize here that they are some important experimental facts in 
favor of the existence of this supersolid state at low temperature in He'', in particular: no NCRI is seen for He^ [55] : 
a (small) thermodynamical signature of a phase transition has been measured for the heat capacity at T ~ O.IK ( [39] 
and see also [40 where it is shown that no clear conclusions can be drawn yet). Amazingly, an increase of the shear 
modulus for solid helium is also observed with a temperature dependence similar to the observed NCRI |41h4^ . 
Notice that this change in the shear modulus cannot explain mechanically the drop in NCRI but could witness the 
occurence of phase transition in solid Helium although it has been observed in solid HcliumS as well! Finally, in what 
could be considered as an "experimentum crucis " Kim and Chan [2^ have shown that the NCRI signal was suppressed 
when changing the topology of the rotating container in order to block a macroscopic superflow. Regarding this 
so-called "blocked annulus experiment" , the only interpretation of the observed NCRI otherwise is by the existence 
of macroscopic superflow, and, to the best of our knowledge, no other interpretation of this crucial experiment has 
been published. 

Given the present experimental context it would be interesting to investigate theoretically the following questions: 

/ - Why does solid helium yield a coherent superflow in a non classical rotational inertia experiment but 
responds as an ordinary elastic solid in pressure/external force driven experiments ? 

II - Why is the superfluid or supersolid fraction so small ? 

III - Why does the superflow density change so widely from sample to sample? and what is the role of the 
He^ impurities in the observed superflow density? 

IV - Why is there an anomaly of heat capacity ? Is this excess of heat capacity related to the superfluid 
fraction ? 

V - Why is there an anomaly of shear modulus ? and how is this anomaly related to the supersolid fraction, 
if it is ? 

The flrst question, /, suggests that, indeed, supersolidity is a rather complex phenomenon: there are two different 
large scale collective motions. A bit like in Landau's two fluid model for superfluid helium, one is the lattice defor- 
mation as in ordinary solids while the second is a superflow, the lattice and the superfluid motion may be realized 
independently. In some situations, depending on the boundary conditions, the system reacts as an ordinary solid, and 
with a superfluid type of motion in other occasions. In this sense, because of rotation the superfluid mode is excited 
or imposed by the boundary conditions, while the solid under gradient pressure does not require superflow to reach 
equilibrium. 

Concerning the second point, //, the zero temperature limit gives a "superfluid fraction" or "supersolid fraction" 
about only 1 %, while in superfluids the superfluid fraction is always 100% at zero temperature. In his seminal 
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paper Leggett indicates that the crystahine structure provides a natural lower value for the superfluid fraction at 
T = OK. However theoretical understanding of this point remains a major, if not the major, challenge. 

Although the points ///, IV and F deserve careful investigations, there is yet no fair understanding of their possible 
connection with the phenomenon of supersolidity and we will leave them for future research. 

The theory of supersolids presents also some outstanding issues. Using an argument already presented by Penrose 
and Onsager [14j, Prokof'ev and Svistunov [44] claimed that a defect-free supersolid cannot exist. Despite the 
presentation by these authors of their claim as a theorem (meaning etymologically a statement that must be respected) 
there has been an ongoing research solving numerically by one way or another the Schrodinger equation to prove or 
disprove the existence of superfluidity in a regular lattice, with widely varying not to say contradictory results. Like 
many others in the field, they assume a fixed lattice without quantum fluctuations induced by zero point phonons 
in the lattice, although those fluctuations are crucial (section II-C) to yield number fluctuations in the system and 
ultimately a long range phase coherence. On the numerical side, based on path integral Montecarlo, there is no 
crystal clear conclusion. For instance, Ceperley and Bernu [45j conclude to the absence of superfluid density in some 
numerical model while Cozorla and Boronat j46j obtain a finite superfluid fraction or NCRI. Clark and Ceperley ||7] 
and Boniensegui et al. |48| predict that a perfect crystal cannot exhibit off diagonal long range order, in contradiction 
with the conclusions of Reatto et al. |49| . It is necessary to be cautious with such numerical estimates: path integral 
Montecarlo estimations are difficult in the limit T — >■ 0. They rely on the winding number that seems to be valid only 
in one (or quasi one) dimensional multiconnected system, known not to be crystal in the ground state. Variational 
approaches and the conclusions drawn depend sensitively on the class of trial functions used. Moreover the analysis 
based on the notion of winding number assumes, although somewhat implicitly, a fixed network, although this network 
experiences quantum fluctuations as well. This type of fluctuations yields a well defined quantum phase to the ground 
state, even if there is one particle per site, an essential remark for the quantum coherence of this state which is 
ultimately responsible of NCRI. Finally, let us mention that the disorder itself has been also invoked to explain 
the NCRI signature observed in experiments, through the concept of quantum glass [SDl EI] and/or the motion of 
dislocation network within the solid [42l HSj |52l [53] . However, to date, neutron diffraction experiments on solid helium 
4 show sharp Bragg peaks typical of a regular lattice structure in the thermodynamical conditions where NCRI is 
observed with no evidence at all of a glass like structure. 

In this review, we present a general description of the theory of supersolids. In particular, one of the aim of this 
paper is to reach a clear and coherent understanding of the macroscopic properties of the supersolid state, so that 
it could be tested in experiments. Although mainly of a phenomenological character, we expect that our predictions 
should be pertinent to solid Helium four at low temperatures (below 0.1 K). Next section (section |ll]) discusses the 
macroscopic properties of a supersolid, based on general arguments and on a Landau type approach. Then, in section 
|III| we show how a quantum coherent phase can be present in ordinary solids at low enough temperature. In particular, 
we explain that the non-zero exchange term between particles in a crystal leads to supersolidity. Finally, before the 



general conclusion, section IV describes a consistent mean-field model of a supersolid that exhibits the main features 



observed in experiments and which satisfies naturally the macroscopic properties drawn in section pl| 



II. MACROSCOPIC DESCRIPTION OF A SUPERSOLID. 



In this section, we show how macroscopic approaches and models can explain the remarkable properties of super- 
solidity. Contrary to microscopic approaches that would describe the dynamics at the atom level, they are based on 
thermodynamic variables such as the mean density, the displacement field and the quantum phase. Indeed, an implicit 
assumption beside is that a kind of quantum coherence is present like in superfluidity. Such hypothesis leans on the 
crucial experiment of the blocked annulus which cannot be explained by an anomaly in the elastic properties of this 
solid. We will first rapidly review the Andreev-Lifshitz two fluids model which was the first macroscopic approach for 
supersolidty. Then, we present how a Lagrangian approach can be deduced from general principles giving a similar 
set of macroscopic equations that in the Andreev-Lifshitz model. Moreover we obtain from this model the consistent 
set of boundary conditions that apply to these macroscopic approaches. It explains in particular the counterintuitive 
fact that rotation induces a superfiow but not a pressure difference. Finally, we sketch the relevant properties and 
some specific solutions of this macroscopic description of a supersolid. 

Let us first emphasize that the present theories are based on the mismatch in the the local number density of 
particles p due to quantum fiuctuations: in ordinary (classical) solids an elastic deformation induces a change of 
density 5p given by the constitutive relation: 

— = -V • u, 
P 
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with u the displacement field. However in a quantum solid, quantum fluctuations give rise to additional elastic 
deformation so that one can expect: 

P 

where gi^" is a quantity (given for the moment) of quantum origin that represents the superfluid or supersolid fraction. 
It is a dimensionless number found to be of the order of 10^^ or 10^'^ in the experiements done on solid Helium. 



Andreev-Lifshitz two fluids model 



In the late 60's, Andreev and Lifshitz [TS], following Landau's approach for the description of superfluidity, developed 
a set of macroscopic equations that would describe the dynamics of a supersolid. These equations provide the 
dynamical evolution of the mass density, the superfluid velocity, the elastic deformation and the entropy transport 
thanks to the total (that is the superfluid plus the ordinary solid part) energy and momentum densities: 



j = rnQ'^'v'' + ■m{p - q'"')u 



(1) 
(2) 



Here £o{p, s,eik) is the internal energy of the body that depends explicitly on the number density p, the entropy 
density s, and the elastic strain eik = {diUk + dkUi)l2. The supersolid density, gi"*, is simply a parameter here. As 
usual, the conservation equations for the total mass, momentum and entropy give |151 154j : 



dtP+^ ■ (j/m) = 
dtji + dkTik = 
(9fS + V • {su) = 0, 



(3) 
(4) 
(5) 



where Tit is the momentum tensor that can be deduced from the energy balance (see below). Additionnally, the 
superfluid velocity satisfies the Landau equation jlOj : 



1 



0. 



Thus the unknown quantities Tik and (f in previous equations are derived from the energy conservation: 

dt£ + V ■Q = Q. 



(6) 



(7) 



After a straightforward calculation, one readily obtains a well defined closed system of partial differential equations, 
using the Einstein's convention of summation for the indices: 



1 s •^2 

m I 

%k = mg''''vfvl+m{p- g'''')uiUk- Sq - Ts ~ fip + ^ g'^^v" - 



Qi = mg" 



P_ 
m 



V ■ u u 

2 



sT+{p-gn{'^u 



+ p] ]ui + aikUk, 



where the temperature is T = the chemical potential is p — and the stress tensor is (Jik — 



the macroscopic set of equations finally reads: 



(8) 
(9) 

(10) 
Therefore 



dtp+V-{g''v' + {p-g'')u)=0, 



0, 



mdtiip - + mg''vld,v'^ + d,{m{p - g'')uiUk) - dk[£o -Ts- p{p - g'")] + dia,k 0, 

dtS + V ■ (sit) = 0. 



(11) 

(12) 

(13) 
(14) 



6 



The first equation ( 11 ) is the usual equation of mass conservation, while equation ( 12 ) is an Euler like equation for 
the superfluid velocity v'^ . In particular, if the initial superflow is not of a vortex type then the superfluid velocity 
displays a potential flow, v'* = with $ a dimensionless field which is eventually linked to the global phase of 

the N body wave function of the system (see below). Equation (13) is the usual equation for the elastic response of 
the system: if no superfluid is present (p'*'' — 0), then one recovers the usual equation for elastic dynamics. Again, 
the last equation ( 14 1 stands for the entropy conservation. Boundary conditions have to be specified for solving this 
set of equations: surprisingly, such conditions are lacking in Andreev and Lifshitz approach and we will deduce them 
from the Lagrangian approach described below. Let us emphasize at this stage that this model opens the way to 
a simple explanation of the apparent paradoxical behavior observed experimentally, that is a nonclassical rotational 
inertia fraction in the limit of small rotation speed but a solid-like elastic response to small stress or external force 



field: the implicit coupling between the superfluid dynamics and the elasticity present both in the equations (11 and 



12 1 might be so that the solid can react very differently depending on the external constraints applied, as we s 
in Sections HTP] and IIV C 41 



lall see 



B. Lagrangian Continuum description of a supersolid 



In this section we present a derivation of the equations of motion for a supersolid in the general framework of 
continuum mechanics and based on general principles. It is based on a Lagrangian formalism, neglecting dissipation. 
The final result includes the equations of elasticity together with the equation for the superfluid component. This 
is valid under the general conditions of applicability of continuum mechanics: in particular the perturbations under 
consideration are long wave (practically much longer than the lattice size) . The form of this Lagrangian is constrained 
by various symmetries, the Galilean invariance and the invariance under rotation and translation. This leaves only 
one possible form, where the microphysics enters only in the values of various coefficients like the Lame coefficients 
of elasticity and the superfluid density p^". The power of this Lagrangian formulation, proposed also by Son |55j . is 
that it does not require any explicit microscopic model but only laws of symmetry. One of the best related example is 
the theory of superfluidity by Landau, free of any explicit connection between the phenomenon itself and its detailed 
microscopic explanation. Besides the general symmetries. Landau made use of the very deep remark that the chemical 
potential is conjugate (in the quantum sense) to the number of particles. A similar property is used below for the 
supersolid, but it has to be changed in order to take into account that the number density of particles may change 
either by compression/dilation of the lattice or by changing the density of the superfluid component. 

This derivation of the equations of motion from a Lagrange principle provides a natural derivation of the boundary 
conditions which are crucial here because they explain why superflow cannot be induced by a pressure gradient in a 
supersolid (contrary to the case of a superfluid) although one is generated by rotation, as explained above. Although 
the model is very much inspired by the deduction by Landau of the equations of motion for a superfluid as a mixture of 
normal fluid and of its superfluid component, it differs from Andreev and Lifshitz approach because of the Lagrangian 
formalism. Here the "normal" fluid is somehow the crystal lattice, although the superfluid is very much like the 
superfluid component of the "mixture" superfluid/normal fluid. Thus, the dynamical equations are obtained through 
the variation of the action S ^ J Ceffdtd^r, where the effective Lagrangian density reads [551 [57] : 



^eff = 



9$ 



fi2 



^ ¥ D^y ^^\P)^ -^Kklm^ik^lm (15) 



In the above expression we have used the so-called material derivative defined naturally as 

Dm du h 



dt 



-V$ • Vu. 



It is the rate of change of a quantity along the flow lines defined by the velocity potential $. This derivative makes 
the rate of change Galilean invariant. As before. 



1 



dxk 



duj 
dxk 



is the strain tensor of the Hookean elasticity theory, Ui being the displacement field. The quantity p|| is the supersolid 
density tensor, which is in general a symmetric matrix. We shall keep it aside for a while, but for more specific 
applications and also for the sake of simplicity, we shall often reduce it to an isotropic form gf^ = g'^'^Sik- Everywhere 
the indices are for Gartesian coordinates. These equations of motion involve the displacement field tti, as in usual 
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solids, the number density p and the quantum phase <J>, related to a ground state wave function. All these three 
quantities are functions of the position r and the time t. 



This expression of the Lagrangian density given in equation ( 15 ) is explained as follows 



1) The part representing the change of energy due to the elastic deformations of the crystal is the same 
as in the standard mechanics of solids. In the limit of small deformations, the final result is minus the 
(positive) quadratic form —\\ikim^ik^im, therefore the elastic stresses are quantities proportional to the 
gradient of the displacement field Ui, with coefficients of proportionality called the Lame coefficients. The 
number of independent Lame coefficients depend on the symmetries of the crystal and is at least two. 

2) From the general principles of quantum mechanics, v"^ = ;^V<i> is the superfluid velocity of the system, 
due to a Galilean boost of the ground state wave function. But in the present case the crystal lattice 
may have a different velocity than the superfluid component. This velocity difference is proportional 
to (V$ — j^%f)- If this difference vanishes (that is if the lattice and the superfluid component move 

at the same speed) the kinetic energy of the motion is just the usual ^pi^^)^ where p is the total 
number density. If only the crystal lattice moves, the velocity ^ is zero and the kinetic energy is just 
half the mass density of the crystal lattice times the square of the lattice speed. Those two possibilities 
(global uniform speed or displacement of the lattice only) are well taken into account by adding the two 
contributions, as done in 



'2m 



p(V$)'-(p5.fe-^.,«fc^) V$ 



m Du\ 



m Du 



This is because the difference between the total density p and the superfluid density tensor, gf^, {pSik — 9tk)^ 
is the mass density of the crystal lattice, a rank two tensor for a crystal in general. This density is multiplied 
by the square of the difference of the superfluid velocity minus the velocity of the lattice in order to make 
this term exactly zero if the two velocities are equal, as it happens when the full system (superfluid part 
and lattice) moves at the same speed: in this case, all the kinetic energy is given by the contribution 
^p(V4>)^ to Ceff so that the other contribution to the kinetic energy must cancel if V<I> — ^j^- 
Notice that this expression of the Lagrangian density can be derived rigorously for the mean field model 
of supersolid [Ml [57] , as introduced in section IV 



3) The macroscopic (or averaged) quantities £, Xikim and gf^ are respectively the internal energy of the 
solid, the Lame elastic constants and, last but not least, the "superfluid" density tensor. These expressions 
do depend, in principle, on thermodynamical variables, crystal structure, etc, and these functions should 
be given by experimental measurements. In particular, we shall explicitly express that the superfluid 
density is a function of the local density number p and it will be often considered further on within an 
isotropic symmetry gf^ = g^^{p)5ik 



Finally, notice that this Lagrangian is Galilean invariant and that the Hamiltonian of the system writes: 

n 



^t^^ + Ut ■ - Ceff, 



Sut 



+ £ + -Kklm^ik^lm- 



(16) 



Taking the variations of £ as a functional of p, $ and u, the final result is a set of coupled partial differential 
equations for the those fields [S51 [57] : 



m 



dt 



dt 

{p^gn[nk--dk<i> 

' TO 



TO 

dxk 



OXk V V TO 



TO 



d 



0. 



(17) 



+ T^{kklmeim) = (18) 

OXk 



(19) 



Here, for the sake of simplicity and further on (except when explicitly stated), we have assumed isotropic macroscopic 
properties for the supersolid tensor gf^ = g'^'^Sik- To these equations the least action or Hamilton principle provides 
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also the boundary conditions, which were lacking for the Andreev-Lifshitz model. In the stationary case for free 
boundary condition it writes: 



Dm. 
Dt 



h 5$ 

m dxi 



dxk 



on dn 
ondfl, 



h Dt 

where is a normal vector to the surface. Similarly in the case of rigid boundary conditions: 



Ui — on dn 
9i$ = on dfl, 



(20) 
(21) 



(22) 
(23) 



Equation ( 18 1 is the usual equation for the elastic response of the system. Equation ( 19 1 is a Bernoulli like equation 
for the superfluid velocity v" = that gives back the usual Landau equation for the superfluid component 

dtv^ = — V^*, ^' being a potential defined directly from (19). Notice that equation (17 1 reduces to the familiar 
equation of mass conservation for potential superflows whenever gf^ — > pSik, corresponding namely to a T = 
superfluid. 

Finally, we can deduce from these boundary conditions the set of boundary conditions for the Andreev-Lifshitz, 
first for free surface (FS) boundary where all fluxes vanish: 



rng'^vtvl 



+ m{p- g''')u,Uk- (So-Ts- pp+^g''''{v' - uf^ S.,k- cr,k 4 = on 

{g""v' + {p- g"')u) -6 = on 
sii ■ e — on 

where e is a normal vector to the surface. On the other hand, for rigid boundaries (RB), it gives 



on 



Ui 

■ e 

s 



on dfl 
on dfl, 
on dfl 



(24) 

(25) 
(26) 



(27) 
(28) 
(29) 



The free surface (FS) boundary condition (24 1 expresses the balance between the normal stress with a kinetic 
(superfluid) pressure. In particular the zero normal stress is recovered when no superfluid phase is present. For 
rigid boundary (RB) condition (27 1 no displacement is allowed at the boundary (obviously, in the case of moving 
boundaries, the condition adapts such that the crystal moves together with the boundaries). The boundary conditions 
for the phase (25 and 28) correspond to a zero flux of matter condition at the boundary (again easily extendable to 
moving boundaries). Conditions (26 and 29) account for the entropy fluxes at the boundaries. 



C. Superfluid Properties: Sound Waves, Non-classical Rotational of Inertia 



Given the general model obtained above, where the superfluid tensor g"'' and the elasticity parameters have to 
be derived from a microscopic model or from experimental measurements, we will present below its characteristics 
properties. In particular, we will explain how the supersolidity appears in the NCRI, we will investigate the long-wave 
perturbations dynamics and the existence of quantum vortices and permanent currents. 



Superfluid density 



The mass conservation equation (17) defines a current reminiscent of Landau's "two-fluid" expression for the mass 
current, at least up to linear order 



Qik^k 



(pSik - gli)uk 



with = — V^. This quantity plays an important role in the NCRI that we shall describe below. Being a 
fundamental quantity for the experimental observations we shall consider in more details how it is possible to compute 
it starting from a microscopic model. 
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2. Non-classical Rotational of Inertia 

As suggested by Leggett [21] , supersolidity should clearly show up via an Andronikashvili kind of experiment of 
non-classical rotational inertia (NCRI), namely a difference between the mass carried under rotational motion and 
the total mass. Experimentally, this is measured with great accuracy by the frequency of a torsion pendulum. Indeed 
let us suppose that the wall of the container of volume f2 rotates with an uniform angular speed uj. For low angular 
speed the crystal moves rigidly with the container ii = uj y. r without any elastic deformation. The densities p and 



g^'^ (where p is the number density) being constant in space, equation (17 1 simplifies into 



V^$ = in n with V$ • e = (m/;i)(a; x r) • e on dVt. (30) 

It corresponds to the classical problem of the flow of an inviscid and incompressible (perfect) fluid inside a rotating 
container. It has a unique solution, and exact solutions can be computed for various geometries in 2D ^58j. For 
instance, as was shown by a student of Kirchhoff, a rotating ellipsoid yields a potential flow that is a linear function 
of the coordinates, a solution given in the book by Lamb 59 1. The effective or measurable moment of inertia comes 



directly from the energy per unit volume of the system (16). In the rotating case it = a; x r the energy (16) can be 
written as 

where Iss is the zz component (the angular velocity ui being along the z-axis) of the inertia moment (we use here an 
isotropic supersolid density tensor): 

Iss = mg^^Ipf + m{p - g'"')Irb 

with 



In 



with $ solution of 

= in n with V$ • e = (z x r) • e on dft. 

This number depends on the geometry only [S51 [50] . The constant Irb is also a geometrical factor corresponding to 
rigid body rotational inertia (x&y orthogonal to the axis of rotation) 



Irb = / {x^ + y^)dr. 



The relative change of the moment of inertia whenever the supersolid phase appears is (here Irb = niplrb) 

(/.. - Irb) g^' (3^^ 



Irb P \ 2^7 



--rb 



Because Ipf < Irb, one has {Iss — Irb) /Irb < as expected and observed experimentally [TJ[51[57]. The NCRI fraction 
disappears when the solid returns from supersolid to the ordinary solid phase {g^'^ — ?> 0). 

3. Sound waves of isotropic lattices 

To characterize the dynamics of this system, we will investigate the small perturbations around a non-deformed 
(m = 0) and steady (V$ = 0) state of average density p using the linearized version of ( |17|18|19D . Before carrying 



the usual small amplitude perturbation analysis, we need to account that an ordinary compression in a solid changes 
the number density of particle in an obvious way, because of the deformation itself, i. e. p ~ p + pV • u + Sp. Finally, 
the linearized system of equations for u, 6^ and Sp reads: 

dt 



m 



_y^s^ - (V • m) =0 (32) 



h—+r{p)Sp = (33) 
m(p - g'')-^^ (dtu - -VS<i>] - (A + A*s)V(V • u) - fisW^u = 0. (34) 
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Here we have assumed an isotropic solid so that the elastic term in (181 simplifies into: 

d / dui 
dxk \ ox„ 

where A = -ftT — is the second Lame coefRcient, and K and /i^ are the compressibility and shear modulus of the 
solid. 



d ( dui \ , . , „ ^2 



Firstly, taking curl of equation ( 34 1 , we obtain that the shear waves are decoupled from other modes and propagate 
following ( n7 = V X u) 

- Q'')-^m - ^.V^tn = 0. 
As fist noticed by Andreev-Lifshitz the shear mode velocity depends on the supersolid density, namely 



^^''---^^(p.g.s)- (35) 

It is thus tempting to deduce from the dependence of the shear mode velocity an alternative measure of ff^^ and 
also to explain qualitatively the increase of the effective shear modulus (extracted experimentally from the shear 
mode velocity) observed at low temperature [35] simply through the growth of the superfluid fraction. However, 
although the increase of the shear mode velocity and the trend of its dependence are correct, such comparison fails 
quantitatively by about one order of magnitude since the shear modulus is observed to vary of more than 10 % for 
solid Helium samples for which the NCRIF is expected to be of the order of one %. Recent results [HI |33] suggest 
that such shear modulus variation could be also due to dislocation motion. However, up to now, let us remark that 
NCRIF and shear modulus have never been measured on the same sample simultaneously. Given the large variations 
of NCRIF between experimental conditions, such simultaneous measurements would be important. 



In addition to this shear mode, taking divergence of equation (34), we obtain for the elastic compressibility (5^ = V-tt 

m(p - [i-M - - V^^*) - (A + 2m,)V2JV = 0. (36) 

at \ at m I 



Together with ( 32 1 and ( 33 ) , it form a coupled set of equations for the compression and for the phase (Bogoliubov-like) 



modes. The dispersion relation for these modes follows from a linear eigenvalue problem for the variables 5p, (5$, and 
dij} in the Fourier space (with S"{p) = j^) : 

E"{p) l.hLu =0. 

{ihuj{p ~ g"')k'^) {{X + 2ps)k^ -m{p- g'')uj^) J \ Si^k J 

Solving this linear system provides a simple algebraic equation for the dispersion relation, from which we can easily 
see that it is linear, so a; = vk, where the wave speed v is given by the roots of the determinant of the matrix, so that: 



2_ K+lfi, (^^ 1^ Ag-^£"ip){p-g^^) 



V — 1 ± /I 

2m(p-£)-) y ]] K + ^p, J 

In the limit p*** — ?► 0, which is the most realistic situations, the two propagation speeds are: 



mp \ p \ K+\pj) 



and 



Introducing the classical longitudinal elastic wave speed ck = \J ^^'^ noting that V2 is the Bogoliubov speed 
for a weakly interacting Bose gas of density g^^ , we remark that: 

S 



Therefore, the second mode, V2 (related to the phase mode), disappears at the supersolid-ordinary solid transition 
while the first mode describes simply the ordinary compression waves. 
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D. Simple steady elastic solutions of the macroscopic equations of a supersolid 



On the other side, we wiU show here that pure elastic behavior can also be observed in such solid under external 
constraint, with no macroscopic quantum phase and thus no superflow. Indeed, this is quite natural since classical 
elasticity of solid is retrieved in the dynamics when no superflow is present. 



For the sake of simplicity we shall assume an isotropic solid so that the elastic tensor term in ( 18 ) reads 



where A 
helium. 



K 



fig is the second Lame coefficient, and K and /Xg are the compressibility and the shear strees of solid 



1. Gravity driven flow 

Let us see if a flow is driven by gravity inside a simple domain as shown in Fig. [T]-a. One has to solve the equations 
( |17|18|19[ ) adding a bulk force ;)er unit mass —gz to ( [18| ) and satisfying the boundary conditions that u = 0. It is 
straightforward to observe that steady solutions are satisfied with an uniform phase, that is a zero superfluid velocity: 



with a steady displacement of the lattice: 



u = 0. 



The equilibrium condition is the elastic equilibrium of a body in gravity for the vertical displacement (u^) in terms 
of the transverse variable (x): 



fi/m + gz 



-mpg 
Ct. 



(37) 
(38) 



where the second condition follows directly (19), assuming that p is constant (or that its variations can be neglected) 
and introducing the chemical potential fi 

Solving the equation (37) for the displacement with fixed boundary conditions, the second equation leads to the 



variation of a chemical potential: fi = fio — mgz. 



x{x — £), 



mpg 
2ps 
/io — mgz 



(39) 
(40) 



Where i is the width of the channel. In conclusion, no superflow is formed in this configuration: the response to 
imposed gravity is a strained lattice without superflow, while, on the other hand, the chemical potential varies with 
the potential energy. 



2. Uniform stress 



We shall focus now on the case of an uniform stress over a supersolid (Fig. [T]-b) . Looking for steady state solutions 
of ( 17|18|19 ) one notices that 



q'^v' + {p- q"')u = ct 
m{p - e'")(ui - 0< - <J,k = Ct 



(41) 
(42) 



Where the first relation (41) is the conservation of mass, and the second one (42) tells us about the conservation of 



momentum and relates in a quite interesting way the elastic stresses with the superfluid velocity. Because the 
symmetries one has that = & it = 0. So that at the end 



fife = Ct. 
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FIG. 1: 



3. Uniform strain 



Similarly, the case of an uniform strain (Fig. [I]-c) presents a no superflow solution because ~ Q 
symmetry. The deformation is simply a uniform compression: 



It = by 



Ae 



4- U-tube 



The U-tube experiment [201 IE2] (Fig. [l}d) combines all the previous cases in a single one: indeed, the left and right 
vertical channels under an imposed gravity while the horizontal tube in between is not submitted to any stress. 

In the vertical tubes, the equilibrium state is given by (37) and (38), while eu 7^ in the horizontal tube and the 
equations of equilibrium read: 



dxk 



Q^'v' + {p- g'')u = 



£\p) = p 



Ct. 



(43) 
(44) 
(45) 



Remark that the chemical potential is constant in the horizontal channel. One can easily see that again a pure "elastic" 
solution with no superflow can be found by combining the solutions obtained above for the three different "bricks" of 
the U-tube. We thus emphasize that the lack of superflow for a pressure or gravity driven U-tube experiment cannot 
prove the relevance of disorder to the observed supersolidity of solid Helium: indeed, no superflow can be induced in 
a good quality crystal by pressure gradient, with or without defects, at least at small stress. 
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III. A POSSIBLE ANALYTICAL APPROACH OF THE REAL DENSE SUPERSOLID 

As already noticed above in section |l] the blocked annulus experiment [2] represents a strong argument in favor 
of the existence of a supersolid state for Helium 4 at low enough temperature. Let us emphasize here that this 
experiment rules out explanation(s) of the observed Non-classical inertia (NCRI) by anomalies in the behavior of the 
torsion pendulum used to put in evidence this phenomenon. It rules out also explanations relying on anomalies in 
the elastic behavior of solid Helium. In superfluid liquid Helium 4 NCRI was observed long ago by Andronikashvili, 
but with major differences with supersolid Helium 4 : 

i) In the liquid, the density of the superfluid component is a function of the thermodynamic parameters: pressure 
and temperature only i[90]. On the contrary in solid Helium 4 the observed density of the superfluid component 
depends on non-thermodynamical parameters, that may be roughly categorized as "crystal defects" because they 
change from one experiment to the other under the same nominal values of the thermodynamic parameters. 

a) In the zero temperature limit the whole liquid becomes superfluid. On the contrary the superfluid component 
of solid Helium 4 remains relatively small at the lowest accessible temperature. 

Hi) The density of the superfluid component, always quite small, tends to zero as the pressure increases in the 
crystal. 

Each of the points made above requires explanation, which we shall attempt to give below. To set the stage, we 
shall outline first some basic ideas about the supersolid state, following the general schema initiated by Onsager |14j . 

A. Phase coherence in a quantum crystal 

We shall introduce first the idea that superfluidity in a crystal is related to the exchange by quantum tunneling 
of atoms from one lattice site to its neighbour. We show that tunneling, along with the quantum fluctuations of the 
phonons, is essential to ensure quantum coherence and so superfluidity, at least in the ground state. 

A standard result of quantum mechanics (a straighforward extension of a classical theorem by Liouville for eigen- 
values and eigenfunctions of real symmetric linear differential operator) is that the wave-function of the ground state 
of a system of Bosons has a constant phase in the absence of rotation. Therefore, it seems obvious that this system, 
whatever is its quantum ground state, liquid or crystal, has the property of off-diagonal long range order (ODLRO), 
understood as equivalent to infinite range correlations of the phase of the wavefunction. Nevertheless, things are not 
that simple because an argument presented below shows that no such order shows-up in a perfect quantum crystal 
(to be defined precisely) if one limits oneself to a particular entry of the density matrix. This raises several issues 
that we address. 

We define first what we mean by perfect quantum crystal, that is a crystal of identical atoms where each atom 
occupies a lattice site with small quantum fluctuations off this site. This smallness of the quantum fluctuations gives 
the idea of a small parameter, ultimately associated to the site-to-site exchange energy and the Debye theory is valid 
at leading order when this exchange energy is small. In a second stage, we consider the property of ODLRO and show 
that it shows-up for some components of the density matrix. 

The idea of studying dense quantum crystals by expansion in a small parameter goes back to a 1960 paper by 
Yang [53]. Yang looked at the ground state of an assembly of identical quantum hard spheres near close packing. In 
this limit, particles have, by definition, little free space. This increases their quantum kinetic energy which diverges 
at close packing with the power law {pc — p)"^, p actual number density, less than but close to pc, density at close 
packing. The inverse square law is explained by the property that the quantum kinetic energy is like the inverse 
square of a distance, and that this length is like the size of the available space for fluctuations of position of the hard 
spheres, a length scale tending to zero like {pc — p) as p tends to pc by inferior values. For arbitrary potentials one 
can try to extend the approach by Yang and to build up a rational theory of dense quantum solids, outside of a 
problematic numerical solution of the Schrodinger equation for a sizable number of particles [51] 

Given the practical closure of direct numerical access to the properties of many-body quantum system with inter- 
actions, one is left with an analytical approach relying on the existence of a small parameter, following the line of 
thinking started by Yang. In the (generic) case of a smooth potential of interaction, this requires flrst to flnd such a 
small relevant dimensionless number, other than the density difference (pc ~ p) of Yang and to expand the relevant 
quantities in the limit where it is small. The zeroth-order theory in this approach is the Debye theory of quantum 
lattices. The quantum Debye theory of solids accounts well for the thermodynamics of crystals, at low temperature 
at least. One expects it to become invalid if the thermal and/or the quantum fluctuations become large enough to 
make the interaction between phonons non negligible. In Debye's theory the zero-point fluctuations are absent. This 
is equivalent to assume that the width of the wave packet of each atom at its site is much smaller than the period of 
the lattice. As well-known such fluctuations preclude the existence of inflnite range positional order in ID. In higher 
dimensions the quantum fluctuations do not destroy the long range order, but make possible the overlap of the wave 
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functions of neighbors on the lattice. This leads to a coherent ground-state wave-function extended all over the lattice 
that is ultimately responsible of ODLRO. Let us put those remarks in an analytic framework. 

The small parameter we are going to introduce is the ratio of the width of the quantum wave packet near each 
site to the lattice size. This small parameter (the (—1/4) power of the quantity A introduced below) tends to zero 
as the number density p increases. As shown below this makes the classical (that is, non-quantum) theory more 
and more exact in this high density limit, and explains that the superfluid density tends to zero when the density 
of the supersolid increases (something hitherto unexplained to the best of our knowledge). When 1/A is small one 
can expand the energy of the ground state explicitly (see equation ( [46| below). We assume a smooth (except at zero 
inter-particle distance) two-body potential. Let us estimate the width of the wave-packet of an atom near one lattice 
site in the field of its neighbors. We do it first for a two-body potential depending with a power law from the distance, 
like V{r) — A and s are positive, the crystal being held by outside pressure. This potential V{t) needs to be 
computed at distances shorter than usual in physical applications, typically fraction of the radius at the minimum of 
a Lennard- Jones potential. Physically, the Lennard- Jones potential increases too much at short distance (like l/r^^) : 
there the dominant interaction is the Coulombian repulsion between the nuclei, far less singular than 1/r^^. This could 



affect significantly estimates of the supersolid density relying on tunneling effects (see IIIB). The present approach 
does not work for the hard-core interaction considered by Yang ^63j. Lastly, the true soft-core interaction at short 
distance could explain why the observed transition temperature to the supersolid state does not change much with 
pressure. 

According to the rules of quantum mechanics the position of each atom fluctuates near the minimum of the 
classical potential at each lattice node. This neglects correlations between fluctuations at different lattice sites, a 
kind of correlation taken into account in Debye's theory through a potential energy term involving cross products of 
displacements at neighboring sites. This does not change fundamentally (at least in space dimensions higher than 1) 
the order of magnitude estimates. At a given lattice site the potential energy due to the neighbors has a minimum and 
varies near this minimum like VsUe ~ T ^ ' "where a is the lattice spacing and q the magnitude of the displacement 
near the minimum. This estimate should include a numerical coefficient of order one depending on the structure of 
the lattice. The first term in this expansion is quadratic in q, the departure from the classical equilibrium position. 

Balancing the potential energy Vsite{q) and the quantum kinetic energy that scales like 2^^, m mass of the He-4 
atom, one finds that the order of magnitude of the zero-point fiuctuations of g is ~ aA~^/'*, where A, a kind of 
dimensionless de Boer parameter, is such that A — ■ In the high density limit and if the potential is steeper 

than quadratic (s > 2, as we shall assume), qq is much smaller than a, the limit we shall consider here. This is 

4 

relevant for most, if not all, real solids, including solid Helium. The dimensionless quantity A^ ~ ^^i~^s'ite(^ — ^) 
is a natural extension of the parameter A for an arbitrary two-body potential. The second derivative Vg-f^{r = a) 
and Young's modulus E of the crystal are linked in such a way that V^lf.^{r ~ a) ^ Finally, the parameter 

= ^^^1- ^ , proportional to the quantity denoted as A for a power law potential, measures the relative magnitude of 
the quantum effects, independently on any assumption on the details of the crystal (two or more bodies interaction, 
etc.). Although \/ Ae is the smallest, about 7.4, for Helium 4, it is still a not too large for solid Hydrogen (12.2) and 
for solid Neon (14.9). 

We shall continue our developments for the case of a power law dependent interaction. The parameter A can be 
used to scale the various contributions to the energy of a dense A^-particle system. Let ip{ri,r2, ■ ■ ■ ,rN) be the 
symmetric A^-particle wave-function. Let us scale rk as afk (tilde being for dimensionless quantities of order 1, k 
particle index), the potential energy as l^(|r"|) = y(a)F(|T'|) and introduce the wavefunction i/)(ri, r2, . . . , rjv) — 
a'^^/^-0(fi, ■r2, . . . , rjv). After the tildes are dropped the energy reads: 

E=^ I I ^ > ;|V.^|^ + A> ■y(\n-r,mAd'^r, 



/ JV N 



where cfi'^r = dridr2 ■ ■ ■ dr^- 

In the large A limit the leading contribution to be minimized in the ground state is the potential energy. If 
one assumes one particle per site, this minimization is exactly equivalent to the one for the classical ground state, 
something done rather easily if one assumes a simple crystal structure and short range forces. The quantum and 
classical problems differ because one does not have to put exactly one particle in each site in the quantum case. The 
wave-function of the full system is a product of normalized one particle wave-functions at a prescribed subset of lattice 
sites among the Ns possible ones. Assuming a proportion p of occupied sites, one adds all contributions obtained 
in this way for all possible combinations of pNg sites. Lastly the wave-function is symetrized under exchange. The 
density of potential energy to be minimized is ^^V{a), the lattice size a being such that p — pa^, p imposed number 
density. The minimization of the energy is done by changing p in the interval < p < 1, p fixed. Decreasing p at 
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constant p obliges to put the particles closer to each other. For a realistic V{r) this increases the potential energy if 
the potential increases rapidly at short distance. Therefore p = 1 is the optimal choice. This defines what we mean 
by "perfect quantum crystal": A should be large and so the optimal configuration is for one atom per lattice sitej92|. 
Notice also that the expansion of any quantity, like the energy for instance, in inverse powers of A is insensitive to the 
quantum statistics (Bose or Fermi) of the atoms. It is at transcendentally small order only that quantum statistics 
become relevant. It is also at an order which is non- algebraic with respect to 1/A that non-classical rotational inertia 
shows up [5S] . 

Let us sketch the principles of an expansion of the ground-state energy in the limit A large. As just said the 
leading order term is the classical potential energy. At the next order, one has to take into account the small 
amplitude fluctuations of the atoms near their equilibrium position. Since the fluctuations are small one can expand 
the potential energy to second order only in the excursion of each atom from its classical equilibrium position. Adding 
now the quantum kinetic energy, one gets a problem of coupled linear oscillators. The diagonalization of the energy 
operator can be done explicitly by Fourier transform and the final result is that the first quantum correction to the 
energy of the ground state is the zero-point energy of the phonons, a quantity of relative order A~^/^, therefore the 
beginning of the expansion of energy per particle reads |64j : 

= Kh^/ma^ (u{p) + A-^/^Eo,ph + ■■■)■ (46) 

The leading order term U{p) is the energy in the classical limit, following a law of corresponding states, and E^^ph 
is the zero-point energy of the phonons, etc. The next order term takes into account the interaction between the 
phonons. This expansion goes on to infinity with integer powers of A~^/^, although logarithms are likely to occur as 
always in expansions generated by quantum theory of interacting fields. However, some physical effects are missing in 
this expansion, most notably the site-to-site tunneling which yields a contribution to Eq transcendentally small with 
respect to the expansion parameter as shown in the next subsection. 



B. Order of magnitude of the tunneling amplitude in the limit A large 



As shown in section [IT] the "superfluid" density g^'' is defined through the superfluid fraction of a rotating quantum 
solid [HI]. This quantity is thus formally obtained through the computation of the variation of the ground state of a 
system of N particles under a small non- uniform change of phase [65] . 

Finally we shall provide a relation between the ground-state wave function of the solid and g^'^ in this large A limit. 
This follows from a chain of arguments already used by Onsager in his derivation of the quantization of circulation in 
quantum fluids [HS]. Consider a real and positive ground state wavefunction: ?/'o(t'i, ^2, ■ • ■ , r^) with energy Eq. We 
shall deal with steady states, so that the energy can be set to zero and one may deal with time independent problems 
only. 

Then one computes the change of kinetic energy by a coherent motion under a small variation of the ground state 
because of a non-uniform change of phase: 

i^{r^,r2, . . . , rjv, i) = e^^^^'-^'+^^'-^^+'-'+^^'-^'VoCri, ra, . . . , rjv) (47) 

where the phase <t>{rj) is taken the same for all particles. Notice that the specific form of the phase of the wave- function 
(47) assumes ab-initio phase coherence, thus superfiuidity. Introducing this Ansatz into the energy ([T]) one gets: 



SS^E^Eo^^ I po(r)(V0)2dr (48) 



2m 

where poif) is the diagonal part of the one particle ground state density matrix 

po(r) EE NRi{r,r) 

with 



Ri{ri,r[)= y dr2dr3...drArV'(ri,r2, ...rAr)V'(ri,r2, . . .Tat), (49) 

where the overline is for the complex conjugation. Notice that the number density poir) > is the true density 
in the system, a non constant periodic function of r if the ground state is a crystal. Because the ground state wave 
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function V'o does not vanish then poir) does not, something true for the ground state of Bosons by a standard result 
of Sturm-Liouville theory. 

According to Ref. [S3] the calculation of the NCRI amounts to estimate the average phase gradient making the 
smallest the change of energy: 5£ = ^ J po(^)(V0)^dr, where V(/) is the gradient of the phase and po{r) is the 
density distribution in the solid obtained by taking the diagonal part of the one-particle density matrix of the ground 
state. We sketch below the proof that, whenever there are wide variations of poir) in the lattice cell (namely when 
the ratio ^/°^p°[^j] is large) , the equivalent superfluid density is the number density at the saddle point in the lattice 
cell. 

The computation of 6£ is done under the condition that the phase (j) changes at large scale at a constant rate, and 
that, given this condition, and a periodic poif), the energy S£ is minimum. For the writing of this energy, it is obvious 
that all the variation of (f) must be concentrated in places where poif) is mimimum. This is even more constraining 
when poi'f') changes widely within a lattice cell. On the other hand, (f>, when it changes from one cell to the next must 
extrapolate between zones of almost constant value to another zone of constant value in the next cell, the variation 
being concentrated in the transition between the two zones of constant phase. Moreover, the contribution to 6£ of 
this transition domain is kept to the minimum by making the transition in the saddle between the two domains where 
the phase is constant because pa{r) is much larger than at the saddle point. According to those considerations, the 
dominant contribution to S£ comes from values of po{r) near saddle points. This can be made more quantitative by 
considering the local expansion of the Euler-Lagrange equation for the optimal phase distribution with a constant 
gradient at large scales. We shall consider now the specific case where the phase gradient is linked to a global rotation. 

Under an uniform rotation the kinetic part of the energy becomes proportional to the square of the angular speed 
of the imposed rotational motion with angular velocity The goal of the calculation is to find the coefficient of o;^ 
in the expansion of the full energy near u; = 0, this coefficient being the observed momentum of rotational inertia. 
Formally this is a standard problem of perturbation. It has been shown that the NCRI can be reduced, for a crystal- 
like ground state, to the calculation of a homogeneized response function. Let po{r) be the space dependent coherent 
density distribution in the solid. For a perfect crystal, this is a periodic and non-constant function of r. The rotation 
of this crystal induces boundary conditions for the slowly varying part of the phase that change the energy. 

Using similar bounds than the one obtained by Leggett [211 IM] one can show that for the given Ansatz used above 
the "supersolid" density satisfy the inequalities 



^ [ min^ po{xX) dr < gl%< [ —, — -r dr 



(50) 



where x is the direction of the boost and C stands for the transversal coordinates. Although these bounds are less 
narrow than those found by Leggett [2T1 [55] both may be estimated in the large A limit via steepest descent. Both 
sides are bounded by the value at the saddle point of po{r) times a constant. Thus a non-zero g'^" is fundamentally 
linked the coherence which comes from the small overlap between the quantum fluctuations at different sites. 

According to the considerations developed above, in the large A limit one expects the density to concentrate in 
narrow domains near the minima of the potential. This yields large density contrasts in the unit cell, because most of 
the cell is filled with low to very low density wave-function. Therefore, as we just explained, the integrals in equation 



(50) is dominated by the value of po{r) at saddle point values in the unit cell. This estimate could be sharpened by 
considering various rigorous bounds of the minimization problem [661 157] -This value of the density at the saddle in 
the unit cell is given by the low density part of the wave-function, off the minimum of the potential, related itself 
to the quasi-classical approximation in the Euclidean case, because this concerns classically forbidden region of the 
trajectories. The general result of this kind of theory is the estimate 
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where 5*0, a Euclidean action ( ie. the action derived from the classical action by changing the kinetic energy from 
/2m to — p^/2m to reach quantum corrections transcendentally small with respect to K), depends on the lattice and 

of the interactions. Notice that the significant quantity in the exponential of the amplitude of quantum tunneling is 

an action, not an energy, contrary to the case of thermally assisted tunneling. 

The action Sq is found by solving the A^-particle Schrodinger equation in the WKB limit for a well-defined 

process of exchange of atoms. We shall show first that the saddle point of the density is on a heteroclinic 

trajectory joining two equilibria (unstable in the Euclidean dynamics). In this WKB limit the wave- function 

ip{ri, . . . , rjv) ^ e~^^o '('"i'---^'"")/''^ where 5Q^''(ri, . . . , r^) is the classical action associated to a trajectory leaving 
a point of equilibrium to reach the point where the action is computed. The long time needed for this exchanges 
yields a cut-off for the frequencies such that the wave- function can be considered as globally coherent. 
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If one takes an arbitrary point in the position space, many trajectories contribute to the amphtude of the wave- 
function there. Along each EucUdean trajectory the action increases with time and the contribution to the density 
decreases with the distance from the equihbrium point. Therefore a general trajectory does not contribute to the 
saddle point density. The dominant contribution is reached when the trajectory is a heteroclinic connection between 
two equilibria exchanging a certain number of atoms between different sites. This corresponds to a stationary action 
and the contributions add to each other to yield a saddle point of the sum of the two sides at the mid-trajectory. This 
saddle point is precisely what we are looking after. The quantity we need to know is the one particle density matrix 
as a function of the position. This is found by squaring the sum of the exponentials and integrating the result over 
all the r2 . . . rpf positions. We consider below a slightly more detailed calculation of Sq in a particular case. 



1. A more detailed calculation in a special case. 

We consider the action associated to an exchange trajectory in the case of a potential V{r) = p-. We shall restrict 
ourselves to a slightly simpler case, the one of two atoms exchanging their equilibrium position at neighboring sites, 
all the other atoms being assumed to be classical points at rest during this exchange, all sites occupied. Let the atoms 
carry indices 1 and 2 at positions ri and r2, the site themselves have index a and /3. The two particles have a large 
repulsive interaction potential AV{ri — r2). The interaction with the rest of the lattice is represented as follows: near 
two neighboring sites a and (3 the potential seen by particle 1 and 2 is a function W{r) = AY^y^a,i3V{r — r^),, with 
all the r-y's fixed for 7 7^ (a, /3) in a first approximation, with two sharp minima at r = and r = rp. This neglects 
possible displacements of the particles other than the one located at and . 

The wave- function for this problem is a solution of the dimensionless eigen-equation (the energy is E = ;^^e): 

eV'(ri, ra) = [V? + V^] 1/; -|- KV{r^ - r^)^; + K{W{r^) + Wir^))^'- (51) 

In the WKB limit, the amplitude of the wave-function under the barrier, that is the tunneling region, is proportional 
to '4>{ri,r2) ^ e~^^o {ri.r^) ^ where S^\ri,r2) is the dimensionless classical action associated to the Euclidean 
dynamical problem with the original kinetic energy — ^ replaced by In this Euclidean dynamics the potential 
wells of the Lagrangian dynamics are replaced by maxima of the external potential. The tunnel factor is found 
by imposing the trajectory to start with particle 1 in well a and particle 2 in well /? and to end up with 2 at a 
and 1 at |93j . In space dimensions higher than 1, such an Euclidean trajectory exists with particles "turning 
around" each other if their interaction potential is repulsive (attractive in the Euclidean dynamics). There is an 
Euclidean trajectory joining states where the two particles interchange their equilibrium positions |94| and in practice 
the pair (ri,a,r2,/3) define the initial condition for the Euclidean trajectory. Fix now a value of (ri,r2), the action 

is S'Q^^(ri,r2) — J ^Vi^q^^ • dri -I- V2S'q^^ • dr2^ = / (pi ■ dri + P2 ■ dr2) where the integral is parametrized by the 
Euclidean time, although pi is the momentum associated to the position r^. Because there are two starting points 

for the Euclidean trajectory, one has to add the two contributions, the overlap being negligible |95j . 

(2) 

Let us sketch the calculation of Sq for a triangular lattice in 2D with a repulsive interaction V{r) = l/r'* in the 
limit s large, quite realistic for a Lennard- Jones interaction increasing like an inverse twelfth power. The lattice size is 
taken as 1 because every other dimensional quantity is absorbed into A j96j . Therefore during the exchange only the 
interaction between the two exchanging particles and with the nearest neighbors is relevant. The trajectories must 
have finite curvature, which imposes that the superposition of all very large forces involved must be such that there 
is no force normal to the trajectories, which would induce a very large curvature. This yields a geometrical equation 
for the trajectory. For a triangular lattice one gets: 

r{e) = (^Vs + sin^e-sine*) 7(273), (52) 

where r{9) is the radial equation with the center at the mid-point between the two equilibria the particles start from. 



those equilibria being located at r — 1/2 and 9 = and tt. The equation (52) is valid for < 9 < n. Using this 



knowledge of the geometrical shape of the trajectory, one finds the analytical expression for the dimensionless action 
^(2) _ \/^3 saddle-point. 



C. OfF-Diagonal long range order or not in a perfect crystal? 



Having shown in the previous subsection the general existence of a non-zero exchange between nearest-neighbours 
sites, we consider now the existence or not of ODLRO in the large A limit for a perfect quantum crystal. 
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Let us consider first with the one-body density matrix and outline the arguments used to "prove" sometime the 
lack of ODLRO in a "perfect crystal". Let '"2, ■■■'T'n) be the full ground state wave function of the N atoms in 



the lattice with one atom per site. The one-body density matrix is given (49 1. The formal existence of ODLRO may 
be understood by saying that, when |ri — r[\ tends to infinity, i?i(ri, r'jjtcnds to the product f{ri)f{r[) with a 
function / ^ 0. This amounts to have that, fixing the position ri, the phase of Ri{ri,r[) does not vary much with 
r[ at whatever distance from ri. Such ODLRO exists for Bose-Einstein condensate, /(ri) being then the mean value 
of the creation (/annihilation) operator at ri. This property of ODLRO, if it is satisfied, is a statistical property, 
and its connection with macroscopic properties like superfluidity is not clear at all: remember for instance that an 
ideal Bose-Einstein condensate (without interaction) is not a superfluid. In this sense, it seems more relevant to deal 
directly with the phase itself, as we shall do further on, and particularly with the relationship between this phase and 
the number of particles |BH1 IM] • 

Coming back to the one body density matrix Ri{ri,r'i), one can argue that it does not show ODLRO {i.e. long 
range correlations between primed and unprimed positions) in a crystal with one atom per site because ri and r[ 
have to lay in the same vacant site of the lattice left once the other sites are occupied by particle number 2, 3, etc so 
that i?i(ri, r[) tends to zero as \ri — r[\ tends to infinity. Recent claims [5¥ on the lack of ODLRO in lattices with 
exactly one atom per site rely on the property of the one-body density matrix we just showed. 

There is an obvious weakness in this argument because it considers only one specific entry of the density matrix 
(the one body), and so cannot exclude that other elements of this matrix show ODLRO, what we are going to show. 
Moreover it considers that the lattice as physically immobile, something incorrect, because it neglects the zero-point 
fluctuations of the phonons. Let us consider first another claim of reference [H] according to which the phase of 
a quantum crystal is not well defined because it fluctuates without bound. Specifically this claim is based upon 
the assumption that, in such a perfect lattice, there is no quantum fluctuation of the number of particles. The 
relation |51 IM] 

5n = 1 

between the quantum fluctuations of the number of particles, Sn, and of the phase 5$ shows that the phase fluctuations 
6^ are unbounded if 6n = 0, implying that there is no ODLRO. This argument seemingly confirm that a perfect 
quantum crystal cannot be a supersolid but it does not apply to real crystals. Actually it omits two types of quantum 
fluctuations of n, some due to the phonons and others to site-to-site tunneling. Even in a quantum crystal with one 
atom per site, the zero-point fluctuations of the phonons, always present, induce quantum fluctuations of the number 
of particles. In ID this kind of fluctuation is so strong that it forbids long range positional order at zero temperature. 
The only meaningful fluctuation in the number of atoms is found by considering what happens in a fixed box of 
volume of order with a size of order L ^ a, a lattice size. This size L must be also much less than the external 
dimension of the crystal, since otherwise the density inside this full crystal does not fluctuate at all. This (obvious) 
remark is not without consequences: practically, to be sensitive to those crucial fluctuations, a numerical simulation 
of whatever kind must be for a system big enough to show finite quantum fluctuations in its subpart. Practically, not 
only N ^ number of atoms, must be large but also N^^^ (see below for the exponent 1/3), which is, practically, much 
harder to reach. 

The number of particles n laying in a finite box drawn on an infinite lattice fiuctuates for two reasons: first because 
of the zero-point motion of the phonons, then, because once tunneling between sites exists (and it does, however 
small it is), there is always an indeterminacy whether a site on the border of the domain is filled or not because 
particles there have a small quantum probability of being on one side of the border or on the other. Consider first 
the contributions of the phonons to 5n, then the one of tunneling. 



1. Number fluctuations due to tfie ptionons. 



Phonons changing n in the box of size L have a wavelength of order L. The effect of the zero-point fluctuations on 
n is found by balancing the zero-point energy of a phonon of wavelength L with the fluctuation of elastic energy due 
to a variation 5n of n in the box of size L (note that the box should be of irregular shape, not a cube or a sphere, 
where the fluctuations in the number of particles depend on the precise way the fixed box is placed with respect to 
the lattice). 

This yields 

In this equation, the left-hand side is the zero-point energy of the phonon of wavelength L with Cg speed of sound 
in the solid. The right-hand side gives the order of magnitude of the elastic energy of this phonon, with E Young's 
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modulus, (1^) strain in the solid, u being the displacement field associated to the phonon. The strain (|^) induces 
a variation of n of order Sn = pjp L^, po mean number density. Once put into the right-hand side of the equation 



(53) this yields the order of magnitude of the mean-square value of zero-point fluctuations of n: 

,hCsP^ 



{SnY 



1 



-paL 



n a 



(54) 



E 

Where = ^^^^ was defined previously. Notice that, in 3D, the fiuctuation of n inside the box are of order 
N^/^ ^ L. The equation (54 1 shows our point, namely that the zero-point fluctuations of the number density do 
not vanish in the perfect quantum crystal. Because they diverge as L tends to infinity, the phase fluctuates less and 
less in this limit. Moreover the scaling (Sn)'^ ~ yields {S(p)'^ ~ fully consistent with the scaling derived by 
balancing an energy proportional to ^ J^g d^r {Wip) ^ L iVip) with the left-hand side of equation (53 1 that is of 
order 1/L at large L |97j . Extending equation (54) to an arbitrary space dimension d yields {Snf' 



the zero exponent is for a logarithmic dependance 5n 
have no well defined position on average, a standard result 



L"^-^. In ID 

log(L) showing that the end atoms of a large fixed segment 



2. Number fluctuations due to tunneling. 

Let consider, as before, a large box of size L drawn on the 3D lattice. The fluctuations by tunneling are due to the 
fluctuations at the sites near the border of this volume. Because those fluctuations are random and not correlated at 
large distances, a fluctuation 6n results from the random addition or subtraction of variables, because there are 
lattice sites near the border. Therefore the magnitude of the fluctuations is like the square root of the number of 
sites near the border, namely like L. For a small probability of tunneling, as expected for real solid Helium 4, this 
small fluctuation is to be multiplied by the (small) square root of the probability of site-to-site tunneling. Notice that 
this source of number fluctuations by tunneling continues to exist in a non fluctuating lattice as well, i.e. in a lattice 
fixed from the outside as the optical lattices for atomic vapors. Therefore, because of those quantum fiuctuations, 
such a system should exhibit long range phase order, at least in its ground state. This could be tested perhaps for an 
atomic vapor with one atom per site of an optical lattice. Therefore, the fluctuation of the number of particles in a 
fixed box of size L is of order L, either because of the phonon zero-point motion or of the site-to-site tunneling. We 
shall use now this estimate to discuss ODLRO. 

The above discussion gives the idea to introduce the quantum amplitude for a state with a given number n of 
particles in a given volume, without specifying which particle is there, contrary to the one-body density matrix 
i?i(ri, r\) that relies on the specification of a particle, something that is not obviously compatible with the quantum 
non-discernability. Let xo.{f) be the function equal to 1 if its argument r is inside a fixed box VL included in the whole 
crystal and zero otherwise, let 5(a,P) be the Kronecker discrete function equal to 1 if a = /3 and zero otherwise, 
a and /3 integers. Introduce now the function Cn('^|Ti, r2, . . . rjv) = '5(X]i=i Xn('''i), "■) that is equal to 1 if there 
are n particles of any index inside a given volume SI and zero otherwise. Thanks to this counting function one can 
derive from the A^-body wave-function 'ip(ri,r2, . . . , t^v) the quantum amplitude (a complex number) associated to n 
particles inside f2 as 

^[n)^ y"d3Ar^Co(n)^(ri,r2...r^). 

The element of the density matrix between states of given number of particles in two volumes SI and f2' is: 

Rn,n'{n,n') ^ j (fi^ rCn{n)Gn'{n')RN{ri,r2 ■ . ■rN;ri,r2 ■ ■ .rN), 

where Rn is the N-body density matrix. The relevance of the fluctuations of n in a given volume J7 for the fluctuations 
of the phase shows-up when considering the Fourier transform of Co(n-), namely the function Dsi('p) such that: 

oo 

Do(¥') = ^Co(n)e™^. (55) 

n=0 

This function Do (1^9) is the quantum amplitude associated to the phase in the volume fi. To get rid of the arbitrariness 
in the phase reference, it is a slightly more convenient to introduce the phase density matrix 

00 00 

Do,a(^,V'') = E E C^,,o(n,n')e™^+™ 

n=0 n'=0 
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Consider now the function Dn^o(</3, f') = Afi{(p — ip'). If n does not fluctuate, An(</3 — f') is just a smooth function 
^moiip-ip ) ^ with no fixed value of n. We have shown that, including in a perfect quantum crystal, n fluctuates in such 

1/2 

a way that ^((5n)^^ is much bigger than 1, therefore Aq{(p — ip') is non-zero in a narrow range of values of ((^ — Lp') 
only. 

This is enough to prove the existence of ODLRO in this system. Take a large volume for f7. From the argument 
just presented, the phase of Co('t-) is a well defined constant. Therefore it does not fluctuate. Take now n ^ n' and £7' 
derived from by a space translation in i?n,n' [n, n') . The lack of ODLRO amounts to the fact i?n,n' (v'j v') tends to 
zero as the distance between O and Vt' tends to infinity. This is true if and only if the phase of Cn (n) and Co' (n) are 
independent, which is clearly not the case because both functions have to have a constant phase and because these 
phases are not independent: they have to be the same phase of a bigger volume encompassing both f7 and f2'. There 
is a slightly less direct argument in favor of ODLRO, that assumes the existence of macroscopic equations of motion of 
the supersohd, as written in section |ll] The phase p of Dj2/((^) is actually the same phase (up to an arbitrary additive 
constant) as the phase entering the macroscopic equations of motion of the supersolid. This is because the phase pi 
is conjugate to the particle number, both as a result of its definition by the Fourier transform as done here or by the 
derivation of the macroscopic dynamical equations by a Lagrange formalism ([S5] and section |ll]-B). Furthermore, the 
equilibrium state of the system, as it follows from the macroscopic dynamical equations, is with an uniform phase, 
the property equivalent to ODLRO in its microscopic formulation. 

Let us come back to the meaning of the statement "one particle per lattice site". Classically it is well defined 
because the counting of particles is unambiguous. Things are different in the quantum case. The number density 
Po{r) in such a crystal is a periodic function of the position and, by definition, the integral of pQ{r) over one cell 
(denoted Vic) , namely J^^ArpQ{r), is the average number of particles in this cell. In general there is no reason 
for this to be exactly an integer in the ground state at a given mean density. The only constraint is that, given the 
number N of particles and the box volume |Otot|, the average number density jfj^ should also be ^ J^^ drpo{r). 
Furthermore the mean number of particles in a box of volume |0| is drpQ{r) = J2^=i ^Rn,n{n, n). The number of 
particles in any box fluctuates because of tunneling and the zero-point fluctuations and, because of that, it does not 
make sense to specify exactly the number of particles per lattice cell. In the ground state at least it does not make 
any sense to specify this number of particles per lattice cell, because, according to the general principles of quantum 
mechanics, this can be done if this number is an observable commuting exactly with the energy operator, something 
that is certainly not correct if exchange between neighboring cells is possible. However this exchange can be small, 
but non zero in the limit A large. In the coming subsection we explore another approach to the same question, namely 
the modelization of a crystal by atoms staying at discrete sites of a fixed lattice. This approach does not take into 
account the vibration of the lattice itself which are essential, as we have just seen, to make the phase of the ground 
state wave function uniform in space. This question will be looked at afterwards. 

D. Representation of a supersolid by a discrete crystal: lattice without phonons 

The theory of solid state makes wide use of lattice models of various kinds, like the Ising model of ferromagnets 
for instance. Of course it is a natural idea to try to model a supersolid by a lattice with "particles" at its nodes. We 
explore below this possibility. 

We consider a model energy operator where particles are constrained to stay on a fixed regular lattice (although 
as we have just seen this assumption of a fixed lattice makes impossible some quantum fluctuations existing in a 
real elastic lattice, see next section for a lattice with phonons). Let / denote the generic lattice index (a set of d 
indices at each node of a lattice of dimension d) . The energy operator representing the energy at leading order is an 
operator £ f{rij) where £ is an energy, large compared to the other energies we are going to introduce, and f{n) 
is a function of the number of particles at the given site. It is a function with a minimum if the site is occupied by 
one particle ( n = 1) and increasing very much if there is more than one particle at the site. It can be for instance 
a function f(n) — n{n — 1)" where k is a large integer. In terms of creation and annihilation operators at the site 
I, n = a^a where and a are the usual Bosonic creation and annihilation operators such that [a^,a] = 1. Indeed 
the energy operator £ fini) commutes exactly with all the nj and its ground state is exactly for the value of 
rij making f{n) the smallest at each site. However this simple property disappears as soon as one adds a tunneling 
piece to the energy operator. Tunneling is described in this simple approach by adding to the energy the operator 
^tunn{o.\aj -\- a^jUi) , the subscript (/J) in being for a sum extended to all pairs of nearest neighbors. 

Define the energy 

n = ^£f{nj) + ^ £tunn{a\aj + a\aj), (56) 
/ (/J) 
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It commutes with the total number of particles, but does not with the number of particles at each site, which cannot be 
considered anymore as specified in this ground sate. Because the total number of particles commutes with the energy 
defined in equation (56), an explicit calculation of this energy requires to impose the number of particles as well. The 
search of the optimal state of this system is slightly complicated. The physically meaningful quantity is the energy 
per particle, given the density of the system. The lattice size is not fixed a priori. It results from the optimization 
of the total energy, including under variation of the lattice size, that is accounted for by changing the energy of the 
single particle energy, under changes of the average density of the crystal. We plan to consider this problem in the 
future. This simple model provides an example of quantum crystal where the number of particle per lattice site can 
be calculated in the (realistic) limit where Sturm ^ S. This can be done by standard algebraic perturbation method, 
considering the tunneling piece as a small perturbation. 

Another remark may be interesting at this stage: simple models like the energy (56) relate at least in order of 
magnitude various physically measurable quantities. It is tempting to do this for the present model. The superfluid 
component is linked to the exchange energy Sturm- This energy yields also an order of magnitude of the transition 
temperature from normal solid to supersolid, this being a bit similar to what London did many years ago to show the 
connection between the Bose-Einstein transition and the observed anomaly in the thermodynamic of liquid Helium 
4 at the lambda point. In the present case, this relates Sturm to the transition temperature. In turn this allow to 
estimate, by a London-like argument, the density of the superfluid component of the supersolid. 



E. Representation of a supersolid by a discrete crystal: lattice with phonons 

The introduction of phonons in this discrete model requires to add one more (quantum) degree of freedom per 
lattice site. This degree of freedom represents the displacement of the atom, if it is there, This section is for the 
purpose of completeness, to show that a lattice model of supersolid may also include the phonons. Suppose we had 
only one atom in an harmonic potential. Its energy is: 

where 5 is a positive stiffness coefficient, and where p and q are conjugate such that [p, q] = ih. The extension of 

this to a lattice of interacting atoms is done as follows. The first term of kinetic energy becomes in the full lattice 
2 

''^i ^ ■ The interaction term should respect the translation invariance of the system: if all atoms move by the 
same amount in the lattice, the energy remains the same. Then a rather straightforward extension of the potential 
energy part of Jis to a full lattice is gJ2(ij) 2 ^J^,/- Because it depends only on the difference between qi and 
gj, it is invariant under a general translation, as requested for the energy linked to the strain in crystal. Therefore 
the full energy operator with the phonon part included reads: 



T-Lphonon = '^{Sf{ni) + -^ni) + ^ i Stunn{a\aj + a]jai) + g 



{qi - qjf 1 /.o\ 
7, '"'I'^'J I (58) 



F. How to reconcile theory and observations ? 



As we have shown, a defect-free crystal of Helium 4 can be a supersolid at low enough temperature if phonons and 
site-to-site tunneling are taken into account. This leads us to discuss how this can be reconciled with the experimental 
findings. 

A rather strong argument for the genuine existence of a supersolid state of bulk Helium 4 is that the observed 
transition temperature from normal to supersolid is fairly independent of the quality of the crystal and depends on 
the pressure only. This rules out any explanation of supersolidity relying on a superflow along a network of dislocations 
and/or grain boundaries, because it would yield a transition temperature depending strongly on the quality of the 
crystal, not what is seen. Another rather strong argument against this explanation by a superflow along a network 
of grain boundaries or in the core of a dislocation networks is that it would need a very high density of defects to 
explain even the small observed superfluid fraction, particularly at the highest pressure where supersolid behavior has 
been observed. We propose to explain the observed sensitivity of the superfluid density with respect to the quality 
of the crystal and to foreign He3 atoms as resulting from the polycrystalline structure of the real crystals. In such a 
structure the superfluid density p^^^ measured by nonclassical rotational inertia (NCRI) experiments, is a function of 
the superfluid density g^'^ in each single crystal and of the way the quantum coherence tunnels from one single crystal 
to its neighbors across the grain boundaries separating the single crystals in the sample. The observed high sensitivity 
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of g^l^ to minute (bulk) concentrations of He3 can be explained by the concentration of He3 in grain boundaries (this 
concentration of impurities in defects like grain boundaries and others is a well-known general phenomenon in real 
crystals). This concentration of He3 atoms likely lowers the grain-to-grain tunneling of Helium 4, including after long 
time scales necessary to reach equilibrium density of He3 atoms in the grain boundaries. The general dependence of 
Sobs with respect to the He3 concentration and annealing time should be rather complex, and not easily describable 
by a single tendency. To check this picture one should measure NCRI in single crystals, to get rid of effects linked to 
the polycrystalline structure of the sample. Somehow this explanation is in agreement with the general fact that it is 
very hard to predict accurately macroscopic elastic parameters of real crystals, because those properties are very much 
dependent upon grain boundaries, etc. The case of superconductors in this respect is quite different of the supersolid: 
in the usual case where the healing length of the superconducting wave functions is much larger than the thickness 
of the grain boundaries, the supraconducing properties are not dependent on the microstructure, polycrystals with 
many grain boundaries or single crystals. 



Below we focus on a model of supersolid valid at T = that is fully explicit and has all the properties requested 
for this state of matter, and that is sufficiently simple to allow in depth calculation of its properties. This is a model 
in the true sense because it cannot be derived from the equations of atomic motion pertinent for solid Helium, it tries 
only to keep the most fundamental properties of this quantum solid to understand its behavior. This model has many 
advantages: the first one is that it is a model of supersolid in the sense of Leggett [21], i.e. it shows NCRI as well as 
an absence of superflow induced by pressure gradient; second, it may be easily implemented in numerical simulations 
where superfluids as well as ordinary solid behavior may be tested; finally, some of predictions may be compared to 
some features observed in experiments. 

Since the work of Kirzhnits and Nepomnyashchii [70^ and of Schneider and Enz |71j in the early seventies the 
transition from liquid to solid Helium has been regarded as a manifestation of an instability as the roton minimum 
in the energy-momentum spectrum touches the zero energy line. Although this idea has been circulating for many 
years, in Ref. [ST] we have proposed a mean field model consisting on the Gross-Pitaevskii equation [711[75][nH] with a 
roton minimum in the dispersion relation, already introduced in 1993 77], that presents a first order phase transition 
to a crystalline state as the roton minimun decreases, but before it touches zero. The difference with the usual Gross- 
Pitaevskii or nonlinear Schrodinger equation (NLS) commonly used for dilute gases where it is a good approximation, 
is in the potential of interaction between atoms that includes now a non trivial dependence on their distance. This 
was the way the original Gross-Pitaevskii equation was written, like: 



where U{\r — r'\) is the two-body interaction potential. At this stage, this model correspond to a superfluid-like 
system at T = OK. Remark that finite temperature can be mimicked for steady situations [TU [7S]. Defining the 
average number density as p = ^ Jy \'ip{r)\'^dr, and noting a the characteristic length of the interaction potential, the 
model has a single dimensionless parameter that governs the type of solutions 



a kind of inverse of the de Boer parameter. This equation is an accurate representation of the atomic reality in the 
mean field limit, namely when the quantum fluctuations inside the range of the potential U{\r — r'|) are relatively 
small. This happens if the average number of particles inside this range is large. This is not the case for real atoms 
where the repulsive forces and attractive forces have about the same range. One could think of electric dipoles to 
enhance the range of the long range part of the interaction. This puts the problem of the ground state in a different 
setting, as the ground state then is made of chains of parallel dipoles, very different of the crystal structure. 
The main properties that may be predicted for this model are: 

i) The ground state of the mean field model is a crystal, that is a periodic pattern (a hexagonal one in 2D and a 
hep in 3D); 

ii) The existence of NCRI; 

in) Four sound modes instead of the usual three in solids; 
iv) Quantized vortex and persistent currents. 

Moreover, in Ref. [SSJ [37] we have developped a theory for the long wave perturbations of the ground state of this 
model of supersolid. We discussed in detail the role of boundary conditions and how to handle steady rotation and 



IV. MEAN-FIELD MODEL OF SUPERSOLIDITY. 




(59) 
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pressure gradient in this model. The mechanical equilibrium was studied under external constrains as steady rotation 
or external stress and the model displays a paradoxical behavior: 

The existence of a non classical rotational inertia in the limit of small rotation speed that do not require 
defects nor vacancies (in full agreement with Leggett's ideas) and no superfiow under small (hut finite) 
stress nor external pressure gradient. The only matter flow for finite stress is due to plasticity being 
facilitated by the eventual presence of defects. 

Our numerical simulations 56J clearly show that, within the same model, nonclassical rotational inertia is observed as 
well a regular elastic response to external stress or forces without any flow of matter, as in experiments (27. .28. 35 , 36]. 

In addition, a new propagating "sound" mode appears besides the usual longitudinal and transverse phonons in 
regular crystals. The speed of this mode is smaller that the usual elastic sound waves speed, it scales as c ~ \/ f"" , with 
jss _ gss I p superfluid fraction ~ 10^**. This slow mode is partly a modulation of the coherent quantum phase, 
like the phonons in Bose-Einstein condensates. Indeed the y/Y^ dependence reminds us the Bogoliubov spectrum. 

The aim of the section below is to study the macroscopical properties of the model given by equation ([59| and to 
compare them with experiments in order to achieve a coherent picture of supersolidity. 



A. The non local Gross-Pitaevskii equation as a model of supersolid 

Our model is based upon the original form of the Gross-Pitaevskii (G-P) equation, which is a non-linear and non- 



local partial differential equation, for the wavefunction of a weakly interacting Bose-Einstein condensate ( 59 1 . This 
is an equation for a complex valued function 'ip{r,t) (this is a complex field, not an operator). 

This equation can be derived from the full Schrodinger equation of many interacting bosons in the mean-field limit 
where the interaction potential is weak and the quantum fluctuations in the range of the potential are small. To have 
stability of long wave fluctuations, the two-body interaction potential, that depends on the relative distance, should 
satisfy that / C/(|r|) dr is finite and positive. Moreover, we shall require also that the Fourier transform 

ik 



Uk = J C/(|r|)e"='-dr (60) 

has to be bounded for all k, and, as we will see later, Uk should become negative in some band- width in the k space. 

The difference with the most commonly used form of the G-P equation is that the potential of interaction between 
atoms includes a non-trivial dependence on the distance although usually this interaction is taken as a 5-function in 



such a way that the cubic term in equation (59 1 becomes simply gip{r)\'il;{r)\'^ . For dilute vapours at low temperature, 
the latter model is a fair approximation of the dynamics, g being proportional to the scattering length for s-waves. 
The authors does not know any experimental situation where this non local Gross-Pitaevskii equation may be applied, 
although it has been suggested recently that a Bose-Einstein condensate with such peculiar interactions could exhibit 
such property [781 179j : whatever the context, such a system would be very interesting for the fundamental physics of 
many-body systems. 

1. Basic properties, invariances and conserved quantities 



The equation ( 59 ) possesses the following properties. 



Suppose V'('''i^)is a solution. Therefore 

-translation invariance : ijj{r + d,t) is a solution as well, with d any constant vector, 
-phase invariance : similarly ip{r,t)e^"' with a any real constant, is a solution. 

-galilean invariance : the "boosted solution" ^{r — vt,t) g;^{™.v-r/h+\mv t/h) ^j^j^ ^ ^j^y constant vector is 
a solution 

Moreover 



Equation ( 59 1 can be put into a Hamiltonian form: 
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where ^p* is the complex conjugate of ip and where the Hamiltonian H[.] is 

Hi^P, ^1 = 1^/ I dr + ^ y y C/(|r - r'|)|^(r)nV;(r')P drdr', (61) 

is positive if U{s) > for all s. 

-The Hamiltonian H, the number of particles N = J \ilj(r)\'^ dr and the linear momentimr P = 
~T / ~ V'^''/'*)da; are conserved by the dynamics with adequate boundary conditions [55] , 

-The Hamiltonian is convex for real values of ip, that is, if one defines E[tp'^] = H[ip,ip] for a real wave- 
function ip, then 20] 

^ ^^2 ^ _ ^)^^2^ < ^^j^2] ^ _ A)i?[^|]. 

After a change to polar variables (sometimes called Madelung transform in this case): ip — -y/pe"^ with 



p and (j) real fields, the Gross-Pitaevskii equation (59) is transformed into the equations of an inviscid 



compressible fluid. They split into the equation for the conservation of matter: 

dp h 



and a Bernoulli-like equation: 



The Hamiltonian takes the form in the polar variables: 

(^|Vp(r-)r + p(r)|V0(r)r) dr + ^ J 



U{\r~r'\)p{r)p{r')drdr', (64) 



Remark. According to the energy (64), the ground-state solution is real (up to a constant phase) and any non 
uniform phase increases the ground state energy. Moreover, in general the ground state cannot vanish. If, however, 
the ground state vanishes at some point as po(r) ^ \r — r*|" with a > one must have a > 2 — D. Indeed a state 



with < a < 2 — D has an infinite (positive) energy, since the energy in (64) diverges as: 



/" ^^|V/o(r)|^drsa finite term -fa^ / |r - rJ^^^dr w finite term + e'(e"+-°"2), 
J 4p(r) Jvir,), 

when e — )■ 0. 

2. Bogoliubov spectrum with rotons 

Given p the mean density defined as the number of particles per unit length, surface or volume in one, two or three 
space dimensions respectively, the homogenous and steady (up to a phase frequency) function ipo = y/pe~^~i^*, where 
Eq = pUq, is a solution of the dynamics, where IJq = J U{\r\)dr (more generally the Fourier transform of [/(•) is 
Uk = /C/(|r|)e''-Mr). 

As discussed in next section, one can show that this solution is stable and makes the ground state for small enough 
n. Indeed, small perturbations around this uniform solution are dispersive waves: 

ip(r,t) — ipQ + [UkC ^ " ' + Vke ^ 'je " 

where fc and uj^ satisfy the Bogoliubov dispersion relation or spectrum |12j 



;..,^,/(— +— pC/.. (65) 
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We shall assume that the potential scales as Uq and possesses a single length scale a. As discussed in [5T], the 
spectrum depends then only on a single dimensionless parameter : A ~ (pa^)^^^^%2 — ^^o that is the product of a de 



Boer kind of parameter (that measures the ratio between the particle interaction and a zero point energy) and of the 
dimensionless parameter pa^ . More precisely we define 



A 



-pUo. 



(66) 



Notice that the existence of a single dimensionless parameter in the problem does not take fully into account the 
complex nature of real solids with, at least, two independent dimensionless numbers pa and ^^^^^^^ — Uo- 
The dimensionless Bogoliubov dispersion relation becomes 



huji, = 



-a)A(fca), with loa{s) = \ + As^ud{s), 



where ud{s) = Ug/a/Uo depends only on the interparticle potential and of space dimension. 
For some analytical results and for the numerics we will choose a soft core interaction [77] 

U{\r-r'\) = UoOia-\r-r'\) 

with 6{-) Heaviside function: 9{s) = 1 if s > and 9{s) = if s < 0. 
The Fourier transform of this special interaction potential reads 



(67) 



(68) 



(69) 



with 



( 2aUo 



41- „3 



1- D; 

2- D; 
5-D. 



and 



ud(s) 



3 / sin s 



COS S 



1-D; 
2~D; 
) i-D; 



where Ji{x) is a Bessel function. Note that Uq is the amplitude of the interaction (with units of energy) while 
Uq = J U{r) d^r, with units of energy x volume {Uq oc a^Uo). 

The spectrum represents well a phonon part together with a roton one (see figure [2| . The long wave fluctuations 
are phonons that propagate with a speed Cg = ^VvV. The short wave regime, the roton part, has the following 
characteristics (see figure [2]): 

if < A < As the spectrum Wfc grows monotonically with k; therefore there is no roton minimum. 
ii) for A = As an inflexion point appears at kg and the energy is given by hujs — Tna^^^' 
in) If As < A < Ac the spectrum exhibits a roton minimum k^a — Sr which is an implicit function of A: 

c3 



A 



(70) 



iv) for A = Ac, the spectrum reaches the axis for fcc as at the edge of the phonon branch in solids. A 
reasonable value for fcc is fcc = 5.45/a (ss 2.1A~^ for Hell). This picture suggest that the existence of a roton 
minimum in Hell is, probably, a reminiscence or a ghost of the solid state as we already suggested |61) . 
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1 2 3 4 5 6 1 ka 

FIG. 2: Various spectrum shapes for different values of A = 12, 23.43, 60 and 90.95 in three space dimensions. 



v) li K > Kc the spectrum becomes pure imaginary in a finite band- width in fc, implying the appearence 
of a linear instability, leading to a periodic pattern of density modulation. 

The parameters considered A^, Ac, A^, Ac, /cs&fcc depend explicitly on the space dimensions. In the following table 
are given the values of the listed above for different space dimensions: 





ID 


2D 


3D 


A. 


9.47... 


15.81... 


23.43... 


kga 


2.89... 


3.33... 


3.72... 




4.92... 


7.26... 


9.82... 


Ac 


21.05... 


46.30... 


90.95 . . . 


kctt 


4.08 . . . 


4.78... 


5.45... 



Finally, we notice that this phonon-roton spectrum has a meaning only as a linear dispersion relation around an 
uniform solution. 



B. Ground-state of the Gross-Pitaevskii model 



With a convenient value of A we obtain a Landau spectrum with rotons. As explained previously, if we increase A 
we expect that there is a critical value Ag < Ac^ < Ac for which the system crystallizes, that is has a ground state 
with a density periodic in space. The increase of A may be realized, for instance, by keeping constant the range a 
and the magnitude Uq and increasing the density n. A density increase might be achieved in a physical system by 
increasing the pressure and/or by cooling. Crystallization due to the roton minimum can be expected near the real 
solid phase since solid Helium exists only at non zero pressure. The transition occurs when the roton minimum is 
near the fc-axis for zero frequency. If we use Landau's notation for rotons: hujk = A + |^(|fc| — fcr)^ for k « kr- In 
our picture A, k^ and are non trivial functions of A. However, A decreases and the roton minimum kr increases, 
as A increases. One should also remark that beside the details of the model, the functions A(A), kr{A) &/Xr(A) are 
known, and ultimately, must be determined experimentally. 

We shall see next via an energy argument that an uniform ground state cannot be stable for any A. As discussed 
above, from ( 64 1 , the phase of the ground state is always uniform in space, even when this state shows modulations 
of the density. A physical consequence is that the ground state has zero momentum P. 

For low A the ground state is a homogeneous solution, a superfluid (without positional order). The uniform ground 
state {ip = ^/p), however, cannot be stable for any A, this fact follows from an argument by Enz and Schneider |71| 
and by Yu. A. Nepomnyashchii and collaborators |70| . Consider a small perturbation around a uniform solution 
p(r) — pq + p[r) and 0(r) — —Eo/ht + (j){r). This allow us to write the Hamiltonian in Fourier (64) in a simple 
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quadratic form : 



Amn 



Uk 



\Pk\' 



Am 





2 " 


4>k 





one sees that if 



Amn 



+ Uk < 0, 



then, the uniform sohition is no more Hnearly stable. Therefore a periodic structure is expected at least as the roton 
minimum reaches the zero frequency axis (that is A > Ac). In |71j the possibility of a linear instability was only 
considered, although the transition is subcritical (first-order) in two and three space dimensions [HI1[7D]. Indeed, by 
decreasing the roton minimum A there is a critical value 1100] such a that is, if A < Ac, then the ground state 
shows a periodic modulation of density in space. 

From now on we shall consider the dimensionless form of the non local Gross-Pitaevskii equation, A being the only 
parameter, defined by (66). In short, the length scale in the problem will be a, the wave function tp being normalized 
by the total density %p ~v^' finally the interaction U{s) will be such a that / [/(It"!) dr 
write the Hamiltonian (61) like H/fl = p£ with Q volume of the system and : 



1. In those units one can 



1 



-\ViP{r)fdr+- f f U{\r~r'\)\ij{r)\^\iP{r')\'^drdr' 
|^(r,i)|2d^r. 



(71) 
(72) 



Here is the total volume of the system in _D-space dimension so that the energy density £ converges because the 
double integral is performed on a compact support (or very localized shape) of the non local interaction U{\r — r'\). 

In the following, we shall estimate modulation of periodic solutions in one, two and three space dimensions based 
on a variational approach. 



1. Weak amplitude periodic modulation in ID 

In one space dimension the minimization of the energy leads to a supercritical (that is continuous second order 
transition) from a homogeneous (liquid phase) solution to a periodic (solid phase) solution. The analysis that we 
follow is the standard perturbation near threshold in the study of pattern formation of lamellar structures |81j . 

If A > Ac a weak amplitude developement with a single wanumber selection is possible, then we consider the wave 
function that is normalized in a period A = 2n/kc according to the normalization condition (72| 



1 



Introducing this trial function into ( |71[ ) one finds the energy per unit length: 

1 2kl\A\^ A , 8\A\^Uk^ ^ 2\A\^U2k^_ 



(73) 



£ = 



1 



2(1 + 2|A|2) 2 \^ ' (1 + 2|A|2)2 (l + 2|A|2)2y 
The minimum of this quantity gives the value for the modulus of the complex amplitude A: 

e + AAUk 



|A|2 = - 



2(fc2+A(C/2fc, -4(7fcJ 



(74) 



The structure displays a periodic modulation with a wave-number kc ~ 4.078 .... Therefore, setting k ~ kc into (74) 
one may calculate this amplitude as a function of A: 



— 8sin(fcc) 
fc3(8-cos(fcc)) 



(A-Ac) w0.011(A-Ac 



28 



2. First order transition in two and three space dimensions 



This ground state can be found by minimizing the energy functional ( 71 1 . Wc sketch the sohition in the case of 
two spatial dimensions. As a trial solution we take ipo the uniform solution plus a modulation, the amplitude of the 
uniform part being set to one by a convenient choice of units : 



with the three complex amplitudes Aj such that \ Aj\ <^ 1, and the vectors kj form an equilateral triangle (fci+A;2+fc3 — 
0) with a magnitude \kj\ = (in this case the critical wave number will be simply the roton minimum, actually 
something realistic in solid Helium w 2-K/kr ~ ZA). If we put iPq{x) in (711 and expand in powers of Aj, since 
l^jl <C 1, we obtain: 



^ = ^ = P^iW'il\^^\'-\ + AlAlAl) + \jZ\Ar 

y i=l i=l 



i<j 



(75) 



withAi=!5A|r). 

The ground state is given by the minimum of the energy, i.e., — 0. The hexagonal crystal solution is {Aj 
Rjc'^'^^) such that each of the three amplitudes is equal to: 



1± ^ 1- 



640 



, J = 1,2,3, 



(76) 



and the phases satisfy (^i + (y92 + ¥'3 = (mod. 27r). The positive sign in (761 gives the stable solution and the negative 
sign gives an unstable solution. The existence of this kind of solution constrained by the condition ji^ < 9/640; or 
better as represented in Fig. [3] : 

WAci(fcr) _ 3 



1,2 



8V10 



has a single (degenerate) root in kr- This gives a critical wave number fc^a = 4.55... and a subcritical threshold for 
the control parameter = 37.36... 




FIG. 3: Various spectrum shapes for different values of A = 37.36, and A = 39 in two space dimensions. The parabola is 
the critical curve ^ J^^ fc^, when it touches the roton spectrum at = 37.36... sub-critical transition towards an hexagonal 
pattern is possible. 
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In physical terms the roton minimum is a more natural parameter than A. In physical units — For 

superfluid liquid Helium (fc^ = 1.95A~^) this condition gives the value piIT] : A^^ = hAK and the corresponding energy 
by particle is about EjN « WT^K. We note that this energy is reduced, probably by hybridization, by three orders 
of magnitude with respect to all of the energies in the theory. A^ is bigger in the three dimensional case, since there 
are more possible stable crystalline configurations, such as hcc or hep. In Ref. [61] we studied the hcc configuration, 
in this case we have six amplitudes and the waves vectors form an octahedron [82], each vector participates in two 

equilateral triangles, which produces a stable configuration. A hcc lattice is stable if A < A^ — which gives 

Aci = 1 .ZK and we have an energy by particle about EjN ft! 5 x \^~^K. We suggest that the supersolid phase occurs 
by increasing the density, therefore an estimate of the critical density is possible, in principle, using an empirical 
law, A(A ~ jo), together with the critical gap A^. Unfortunately we do not have adequate experimental data (A as 
a function of n). However, from |83j we can fit A as a function of the pressure (p). The fit gives A = A(0)(1 —p/po), 
with A(0) = ^.7AK and po = 1576ar. We find a critical pressure p^ = 26bar for the transition to a bcc lattice, in good 
agreement with the experimental value. 

Numerically, the ground state can be reached by considering the dissipative version of the G-P equation, called 
the Ginzburg-Landau (G-L) equation that can be understood as a imaginary time evolution of the G-P equation. It 
writes in dimensionless form: 



^ = - V^V - A^(r) J U{\r - r'|)|^(r')pdr' + ^V, (77) 

where is the chemical potential to impose the mean density (or total mass) of the system. The G-P and G-L 
equations have the same stationary solutions and ground states and the dissipative G-L dynamics converges to a local 
minimum of the free energy that can be deduced from the Hamiltonian. Thus, starting with a noisy initial conditions 
and running numerically the G-L dynamics, one can achieve a ground state. We use two types of numerical methods, 
a pseudo-spectral one when periodic boundary conditions can be considered and a finite difference one based on a 
Crank-Nicholson scheme otherwise. Figure Q shows the ground state in 1 and 2 dimensions with periodic boundary 
conditions. A regular ID crystal is observed and in 2D an hexagonal pattern is formed, minimizing the free energy of 
the system. 




FIG. 4: Ground state of the G-P equation obtained numerically using the dissipative G-L equation a) in ID, where a periodic 
crystal of density peaks is formed (here A = 43.2); b) in 2D, it leads to a periodic hexagonal pattern (A = 107). 



In three space dimensions the most stable configuration is the hep crystalline structure as it can be seen on figure 
([5]). In this case the critical wave vectors belong to a tetrahedron, whereas in the 6cc structure each vector belongs 
two resonant equilateral triangle. 
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FIG. 5: Ground state of the G-P equation obtained numerically using the dissipative G-L equation in 3D, exhibiting a regular 
hep crystal (A = 429). 



3. Ground-state in the large A limit. 



The crystal structures with a weak modulation of density have been described in 1-D for moderate values of A, say 
A > . In fact, the ground state defined as the minimum of the energy functional ( 71 ) exists and is unique, because 
of the convexity, for any A. For larger values of A the weakly nonlinear approach cannot be used anymore to derive 
the amplitude of the density modulation. 



For large A, the potential energy in ( 71 1 requires small tp while the mass normalization ( 72 ) forbids ip to be small 
everywhere. Therefore the energy minimization leads to a periodic structure with zones where ^ w balancing zones 
where V' 3> 1. Recently it has been proven by Aftalion et al. [671 a remarkable, in the opinion of the authors, theorem: 
in the limit A — > oo, the minimization of (71) with the restriction (72) is equivalent to a close packing arrangement 
of rods, disk and spheres in one, two and three space dimensions respectively. 

In the particular case of 1-D, it is shown, in Ref. |67j, that in the large A limit ip only in a small zone : 
X e (— <^) of the whole period: (—A/2, A/2), where the Euler-Lagrange condition deduced from (71 1 together with 
(72) leads to the Hemholtz equation in the domain (—S,S) : —ip"{x) = fiip. Finally the minimization of the energy 
gives S and the wave number A of the periodic structure. 

Following this approach, with Sepiilveda ^5] we estimated such a ground state for A ^ 1 (the extension to higher 
dimensions seems natural but the computations are harder). We sketch below the corresponding results. We consider 
the energy and normalization ( 71|72 ) in a single period with the trial function in the unit cell |102] 







. 25 ) 







X e [-A/2,-,5] 
X e [-(5, S] 
X e [S, A/2] 



(78) 



The integrals maybe computed more or less easily because of the interaction term. The minimization yields a relation 
between S, the wavelength A and the dimensionless A.: 



l + 2(5-A = 2 



(A(A-1)2A)^ 



(79) 
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This variational result is in complete agreement with the numerics, a more complete discussion on this problem may 
be found in Reference [86] , 

4- The question on the commensurability in the Gross- Pitaevskii model. 

The question of the number of "atoms" N versus the number of sites Ns and how are they related, if they are 
a relation, needs some discussion in the context of the present model. The problem if the solid is commensurable 
or incommesurable, that is if Ng/N = 1 or Ns/N ^ 1 is, in the context of this model, rather different from a 
commensurability problem in a classical solid. More precisely, let us consider a initial wavefunction with initial 
number of atoms (or mass) / IV'pdr = A'^, being an integer in a periodic box of size L, L x L, L x L x L in one, two 
and three space dimensions respectively. The number density is, therefore, N/L^ . 

In the present model the crystal appears by a spontaneous breaking of the symmetry under translations of the 
uniform state. Neither the sites nor the peaks (or "atoms" ?) were there before. Moreover, as in every spontaneous 
symmetry breaking, the appearance of the structure is a dynamical process governed by transient structures full of 
defects, vacancies, interstitials, grain boundaries, etc. In infinite length domain, it is expected that these defects, 
vacancies, interstitials, grain boundaries, etc. move away the system leading to a defect-free crystal. This process 
may take, however, a very long time. On the other hand, for finite length system, which is the particular case of 
numerical simulations, the formation of a periodic structure needs to match the boundary conditions that makes 
difficult the realization of a free defect system. 

As an example, the "pattern formation" process in one space dimension is such a that the wavelength cannot 
be selected via the most unstable mode Xc — because the wavelength should match the boundary conditions, 
something that becomes more or less irrelevant as the system size grows L ^ Ac. The number of sites is related to 
the selected wavelength via XNg = L. So in practice the ratio of the number of sites to the number of "atoms" (the 
initial normalization) N are not the unity N/Ns = NX/L = pX except for very particular cases. In higher dimensions 
the mismatch of the crystal structure-boundary conditions becomes more complex. 



One of us (SR) has performed 2D numerical simulations of equation (59) in a 64 x 64 periodic box. The interaction 



range was chosen to be a = 8, the mesh size being dx = 1. The normalization condition is for J \ip\'^dr = iV = 36 



"atoms" . Although the simulations of the original Gross-Pit aevskii equation ( 59 1 are formally reversible, an irreversible 
behavior appears naturally that takes the energy from long scales to the smallest ones ~ dx driving the system finally 
to a macroscopic state that minimizes the energy |85j . It is important that the relevant structures are much larger 
than this size. Indeed in the present case the roton minimum is placed at kr ~ 4.55/8 = 0.57 which is much smaller 
than the wave number associated to the mesh size: tt. 

One sees that despite the same initial condition, for the four cases presented in the figure |6j the system reaches 
different commensurability fraction N/Ng both larger and smaller than the unit depending on A. In fact, for larger 
A (that is larger compressibility because A ~ p) the system makes in average more sites of lower mass each. And 
conversely for small A (small compressibility) the system makes a smaller number of sites. We can conclude that 
in the present model the commensurability is not a relevant parameter, as suggested by the mean field origin of the 
model. In fact, despite the different commensurability, figures a), b) or c) and d) are qualitatively the same. 

This question of commensurability could be investigated experimentally. If there is in average a non integer number 
of atoms per lattice site, by neutron (or X ray) scattering, one could see the mismatch between one and the actual 
number of atoms per site. This is surely not an easy experiment, because the mismatch is likely small and it requires 
an accurate measurement of the mass of the Helium crystal unit for diffraction. Perhaps a differential experiment 
could be done by changing the volume under applied pressure. 

C. A model combining elastic and superfluid properties 

1. N on- Classtcal Rotational of Inertia 

As a built-in superfluid model, it is expected that the Gross-Pitaevskii equation exhibits superfluid properties in 
presence of a crystal lattice. This has been shown first in [ST] where the superflow of such supersolid model was 
observed around a cylindrical obstacle with a critical velocity above which vortices where nucleated, similarly than 
for the superfluid model [84] . 

In fact, it is possible to use this model to mimic the torsional oscillator experiment where NCRI has been exper- 
imentally observed [T]. In [S5], we have observed the NCRI effect using a 2D numerical simulation where a square 
sample is put under rotation. Measuring the rotational inertia, we have demonstrated that a part of the total mass 
was decoupled from the rotational motion, as shown on figure ^)). We use a relaxation algorithm (based on the 
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FIG. 6: Various density plot for long-time (1000 time units) numerical simulation of (59 1 with a = 8 in a 64 system size with 
periodic boundary conditions. In all simulations the initial conditions are identical and A'^ = 36. The values of the interaction 
parameter Uo are: a) Uo = 0.5 and A = 56.55 ; h) Uo = 0.7 and A = 79.17; c) Uo = 0.75 and A = 84.82 and d) Uo = 0.85 and 
A = 96.13. The number of sites are, respectively, a) Ns = 34; b) Ns = 36; c) Ns = 36 and d) Na = 40. 



Ginzburg-Landau equation which consist on a imaginary time evolution of the G-P equation as explained above) to 
converge towards an equilibrium state of the system (close to ground state). The numerical simulation are performed 
using a pseudo-spectral method and the system is put under rotation by considering the equation in a constant rotat- 
ing frame. Figure (jjjD)) shows the evolution of the NCRIF in the limit of zero angular valocity as function of A (here 
by varying pUo), indicating that the superfluid fraction decreases as the pressure increases in agreement with Leggett's 
argument that the superfluid fraction should decrease with the minimal value of the wave function j211 [571 [SB] . 

A more accurate way to compute the NCRI can be obtained by considering a small fraction of the sample only (see 
figure ([s])). In such case, neglecting the centrifugal acceleration term (which would induced a correction in on the 
solution of the wave- function) , the rotation can be approximated by a Galilean boost in the azimuthal direction. The 
limit case of a Galilean boost of a ID system can be understood in this context as the rotation of a very thin torus 
of the crystal. With N. Sepulveda, we have numerically computed the NCRI in 1 and 2 spatial dimension using this 
method, allowing an accurate measure of the superfluid fraction in this model of supersolid [551 157] . 

Figure ^a.) shows the NCRI as function of the Galilean boost velocity in 2D, showing a decrease of the NCRI as 
the velocity increases. Moreover, in this configuration rapid variations of the NCRI appear regularly and can be. As 
shown on figure (|9]-b), these sudden decrease of the superfluid fraction is related to the nucleation of quantum vortices 
which are known to lead to the loss of superfluidity in superfluids [55]. 
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FIG. 7; 2D numerical simulation of the dimensionless equation with 128 x 128 modes in a square cell of 96 x 96 units for 
different values of pUo; the potential range is a = 4.3. (a) The NCRIF = 1 — L^cj)/ (Irt) vs. the local maximum speed 
Vmax = u}L/V2 for pUo = 0.069, 0.084, 0.099 & 0.114 (res. A = 74.1, 90.2, 106.3, & 122.4). Here (7^) is the converging inertia 
moment computed numerically for large nUo at w = 0. Note that the jump in NCRIF for pUo = 0.069 (A — 74.1) corresponds 
to the nucleation of a vortex in the system, b) NCRIF at oj = as a function of pUo- We have verified that (a) and ( b) are 
almost independant of the box size. This figure is similar than that in 




2. Quantized vortices and permanent currents 

In general, the existence of a macroscopic quantum phase suggests that quantized vortices and permanent currents 
should be present. Such supersolid vortices" consists of a ±27r jump of the quantum phase around the so-called 
vortex core. It corresponds thus to topological defects that cannot be removed by infinitesimal perturbations of the 
system. Actually, quantum vortices can be observed within a mean-field approach of a quantum solid proposed here, 
as illustrated in figure [TOl Using an energy minimization argument, one shows that the core of the vortex is located 
at a minimum of density of the crystal lattice. 

How could these vortices be observed in real supersolids? Imagine a number of superfluid vortices in a Hell 
superfluid created by a rotation at very low temperature, (typically below O.IK). Then increase the pressure in order 
to solidify Helium. If it is a normal solid phase, no quantum phase is anymore present and the vortices disappear. 
But, if it is a supersolid, then the evolution of the vortices during the crystallization process should be as follows: 

1) The vortices in the fluid are pushed by interaction with the crystal boundary and driven to the container 
boundary; or, 



2) The crystal grows around the vortex to form a "supersolid vortex" (numerical simulations of (59) seem to confirm 
this second possibility). 
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FIG. 9: (a) Supersolid fraction as a function of the boost velocity when the crystal is submitted to a Galilean boost. The 
continuous line shows the component of the supersolid tensor parallel to the velocity while the dashed line is for the orthogonal 
direction, b) Phase of the wavefunction for the boost velocity Vh = 0.4. It shows twice a 2tt phase jump (the phase is defined 
with a multiplier of 27r) indicating that vortices are present in the sample (see figure (|8|). 
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FIG. 10: Numerical finding of an energy minimum with a topological vortex in the middle of the system in the frame of 
the energy (61 1 with a constant total number of particles A'^ = 121, the range of interaction is a = 8 and the size system is 
112 X 112 with Neumann boundary conditions, that is tp{x,y) = if {x,y) are outside the domain. The interaction parameter 
is Uo = 0.75, that is A = 82.94. The right-hand-side figure a) represents the density of the crystal, that is \ip\'^ with Ns = 110 
sites, while the left-hand-side b) represents the phase of the complex wave function tp. A phase jump is clearly visible on the 
left side of the vortex which is located at the center of the system at the end of phase jump. It is interesting to note that 
the phase maybe split in two parts, a slowly varying part that change over the system size and a periodic perturbation with a 
well defined relief of the crystalline pattern as is guessed in the homogenization theory. See section |IV] for more details on the 
model. 



Finally to see experimentally the existence of supersolid vortices, we again melt the supersolid. The vortex cannot 
disappear and thus we obtain a vortex in liquid Hell which could be detected as done in superfluids. Therefore, if 
one could show experimentally that a vortex survives the freezing-melting process, we would have a clear proof of the 
existence of a supersolid phase! 

However, the first possibility might happen and no vortex could be measured while the supersolid phase still exist! 
A way to avoid this problem could be by freezing a circulating superfluid in a multi-connected domain (103) . Similar 
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Gedanken experiments could be considered using permanent currents. A permanent current may be obtained in a 
supersolid filling a multi-connected domain. Indeed, imposing a phase jump of 27r as one turns along a non collapsible 
closed curve in a multi-connected domain we assure the existence of a non-uniform phase as in the case of a vortex. 
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FIG. 11: Numerical finding of a energy minimum with a permanent current in the frame of the non local Gross-Pitaevskii 



energy (61 1 with a constant total number of particles A'^ = 96 the size system is a multi-connected domain composed of a 
112 X 112 square box with a circular hole of a radii of 15 units in the middle, and we use Neumann boundary conditions, that 
is 'ip{x,y) = if {x,y) outside the domain. The range of interaction is a = 8 and the interaction strength is Uo — 0.75, that 
is A = 82.94. The right-hand-side figure a) represents the density of the crystal that is \ip\^ with Ng = 106 sites, while the 
left-hand-side b) represents the phase of the wave function. A phase jump is clearly visible on the left side of the hole. 



Vortices and permanent currents are a non-ambiguous property that would confirm the existence of a supersolid 
coherent state of solid Helium. To date there is no direct experimental evidence of such properties in solid Helium at 
conditions where NCRI exists. 



3. Elastic properties 



The crystal structure that forms naturally in this model exhibit also classical elastic properties. Obviously, the elastic 
response of the system is due to the rapid density variation of the mean field in the crystal. Using a homogenization 
technique ^9] that separates the small scale short time variation at the level of a density peak with the long wavelength 
slow macroscopic modes, we have been able to deduce from this G-P equation the macroscopic model described in 
section |ll]-B |56l 15 7j . In addition, this technique provides an explicit protocol to compute the superfluid tensor as it 
has been checked in details in [87] . 



4- The U-tube experiments visited in the numerics 

In Ref. [5G\ we studied a gravity driven supersolid flow. As early suggested by Andreev et al. |18j an experiment of 
an obstacle pulled by gravity in solid helium could be a proof of supersolidity. Different versions of this experiment 
failed to show any motion [62], therefore a natural question arises: How we can reconcile the NCRI experiment by 
Kim and Chan and the absence of pressure or gravity driven flows? 

As already said in section 



ways under a small externa 



II D our supersolid model (and it seems that supersolid helium too) reacts in different 



constrain such as stress, bulk force or rotation in order to satisfles the equation of 



motion and the boundary conditions. For instance, if gravity (or pressure gradient) is added then the pressure £'{p) 



balances the external "hydrostatic" pressure mgz in equation ( 19 1 while the elastic behavior of the solid of equation 



(18) balances the external uniform force per unit volume mpg. No V$ nor ii are needed to satisfy the mechanical 
equilibria. Moreover, a flow is possible only if the stresses are large enough to display a plastic flow as it happens in 
ordinary solids. In [61] we showed that a flow around an obstacle is possible only if defects are created in the crystal, 
in this sense we did observe a plastic flow, however in the same model we observe a "superfluid-like" behaviour under 
rotation without defects in the crystal structure. Indeed for a small angular rotation the elastic deformations come 
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to order llP' while V$ or u are of order w, the equations of motion together with the boundary conditions leads to a 
NCRIF different from zero. 

We have realized a numerical simulation to test the possibility of a per man ent gravity flow for different values of 
the dimensionless gravity Q — . Let us consider an U-tube as in Figjl2[ The system is prepared for 500 time 

units into a good quality (but not perfect) crystalline state. A vertical gravity of magnitude Q is switched-on and the 
system evolves for 500 time units more up to a new equilibrium state (see Fig. [l2]-a). The gravity is then tilted (with 
the same magnitude) at a given angle. A mass flow is observed at the begining from one reservoir into the other, 
but both vessels do not reach the same level eventually (see Fig. 12 -6). There is some dependence of the transferred 
mass on Q till Q « 0.0005 and the mass transfer becomes negligible from fluctuations for Q < 0.00025 indicating 
the existence of a yield-stress. The flow is allowed by dislocations and grain boundaries and it is a precursor of a 
microscopic plastic flow as in ordinary solids {e.g. ice) and as it is probably observed in Ref. [62] . A microscopic 
yield-stress could be defined by the smallest gravity Q such that no dislocations, defects nor grain boundaries appear. 
In the present model this is for Q < 10~^. 







500 1000 1500 2000 2500 3000 t 



FIG. 12: We plot the three snapshots of density modulations \ip\'^ (the dark points means a large mass concentration) of 
a numerical Simulation of eqn. (591 with Dirichlet boundary conditions with the shape of an u-tube as in the figure. We 
use a Crank-Nicholson scheme that conserves the total energy and mass. The mesh size is da; = 1, the nonlocal interaction 
parameters are chosen as Uo ~ 0.01 and a = 8 (physical constants h and m are 1), finally the initial condition is an uniform 
solution tp = 1 plus small fluctuations, a) A crystalline state is obtained after 500 units of time; b) then a vertical gravity of 
magnitude Q = 0.01 is switched-on, and the system evolves for 500 time unites up to b: the solid makes a neat interface; c) 
Finally, the gravity orientation is tilted in 45°. After 2000 time units the system evolves to a stationary situation c showing 
that the mass flow is only a transient. Moreover, in d) it is observed, after 1000 time units, a stationary situation showing that 
the mass flow is only a transient. 



In conclusion, we have shown a fully explicit model of supersolid that display either solid-like behavior or superflow 
depending on the external constrain and on the boundary conditions on the reservoir wall. Our numerical simulations 
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clearly show that, within the same model a nonclassical rotational inertia is observed as well a regular elastic response 
to external stress or forces without any flow of matter similarly than in the macroscopic model and in agreement with 
the experiments [H . 



Sound as an alternative test for superso 



As explained in [SHIEZIj this model, beside the NCRI property, has an intrinsic elastic behavior and follows the 
macroscopic equations discussed in section [llj-B. In particular, it has been shown numerically there that the quantum 
solid of the model can present a pure elastic response to an imposed external strain due to gravity. Let us also recall 
that such solid exhibits different perturbation modes (sound waves) , a pure shear mode on one hand, and compression- 
phase coupled modes on the other hand. In particular, in the limit of small superfluid fraction, a low-velocity mode 
exist which has the same characteristics than the Bogoliubov mode for superfluid or Bose-Einstein condensate. In this 
limit, the velocity of this phase mode is proportional to the square-root of the superfluid fraction of the supersolid 
(ti| = ^ — ^— ^). The existence of such mode suggests an alternative way of determining the superfluid fraction in a 
supersolid (and by consequence an alternative test for showing the existence of the supersolid state), by looking to 
the resonant frequencies of a "supersolid" cavity. In this paragraph, we will show numerically in the present mean 
field model that the lower frequency eigen-mode of a cavity provides a correct (and thus alternative) measure of the 
superfluid fraction. 

In order to measure numerically the cavity modes we are interested in, we create a low frequency oscillation of 
the system by starting the numerical simulation of eq. ( |59| in a closed cavity of size LxL with a small initial mass 
difference between the left side (0 < a; < L/2) and the right side {L/2 < x < L). More precisely, we use as initial 
condition: 



I - e + noise < a; < L/2 
I + e + noise L/2 < x < L 



with typically e = 0.1 and a zero mean noise of amplitude 0.05. We will vary the dimensionless parameter A by 
varying in the numerics the value of Uq for a fixed. By measuring the time evolution of the mass in the left side of the 
cavity {Ni{t) = dy dx\il){x^y)\'^), we can extract the vibration modes of the cavity using Fourier transform. 



Figure ( 14 ) shows such frequency spectrum and one can identify the eigen-modes of the cavity as the different peaks 



of the spectrum. 



'^'^ ^ o 




FIG. 13: Snapshot of the density modulations a), and the real part of the wave function Reip, b) of numerical Simulation 
of eqn. ( 59 1 with in an annular domain with Dirichlet boundary conditions in the horizontal axis and periodic boundary 



conditions on the vertical axis. We use a Crank-Nicholson scheme that conserves the total energy and mass. The mesh size 
is dx = 1, the nonlocal interaction parameters are chosen as Uo — 0.005 and a = 8 (physical constants h and m are 1), both 
snapshots are at the same instant of time t = 125, 000 time units. Though the modulus seems more or less uniform in the large 
scale one sees the long wave oscillations of the phase of the wave function in the figure 6. 



Seeking to identify the low velocity mode with the low frequency cavity, we exhibit the two smallest frequency peaks 
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as we vary A, by changing Uq and the initial density, and for different boundary conditions (usual conditions is no split 
at the boundary) in table |l] We also compute for these values of the parameters the superfluid fraction bounds using 
the Leggett theory [5T]. Leggett deduces from general argument an upper and lower bound for the superfluid density 
which is easy to compute when the wave function is known (see [55', 'ST for a discussion on Leggett 's formulae). As 
it can be seen in the table, the lower and the upper superfluid fractions are very close one from each other and we 
will take the superfluid fraction as the average of both bounds = (/!" + /^*)/2 ± {{f^^ — ft^)/2. 
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FIG. 14: Fourier transform of the signal Ni{t). 



Fig. [Ts] shows these two lowest frequency modes as function of this We observe that the flrst mode satisfles 
a Bogoliubov-like dependance on the superfluid fraction as expected from the macroscopic theory while the second 
mode does not depend on the superfluid fraction and should be associated to a pure elastic mode. 
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b) 



FIG. 15: The first resonance /i (a) and second resonance (b) versus the supersoHd fraction measured thanks to Leggett's 
iformulas applied to the model as explained in the text. In a) The fine represents the curve 0.032v7^ for guideline to the 
eyes. The two points represented by a ■ are low values Uo a situation where the ground state is homogeneous in space and 
represents better a liquid than a solid. 



From this preliminary numerical study, we suggest that an alternative (and complementary) way to measure (and 
somehow to exhibit) supersolidity could be by identifying pressure and/or temperature dependence of the frequency 
of the lowest eigen-mode of vibration of solid Helium. 
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Uo 
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/i/lO 


±A/i X 10"^ 


f2/W 


f- 


rss 

J + 


p=l 




0.00300 ■ 


38.6039 


1/39 


6.58 


1/16 


0.984752 


0.9983 


0.00325 ■ 


41.8209 


1/41 


5.95 


1/17 


0.913244 


0.920613 


0.003375 A 


43.4294 


1/42 


5.67 


1/17.7 


0.559121 


0.607666 


0.00350 ▲ 


45.0379 


1/43.5 


2.64 


1/18.5 


0.481871 


0.531268 


0.003625 A 


46.6464 


1/44.5 


5.17 


1/14.7?? 


0.422077 


0.490596 


0.00375 A 


48.2549 


1/47 


4.53 


1/15 


0.376049 


0.432442 


0.003875 A 


49.8634 


1/61 


3.97 


1/17.5 


0.33623 


0.408770 


0.00400 A 


51.4719 


1/57 


3.08 


1/19.5 


0.303137 


0.365992 


0.00425 A 


54.6888 


1/62 


2.69 


1/18 


0.247761 


0.311837 


0.00450 A 


57.9058 


1/65 


2.36 


1/19-1/15 


0.203728 


0.269468 


0.005 A 


64.3398 


1/74 


1.82 


1/18.5 


0.141113 


0.206708 


0.006 A 


77.2078 


1/89 


1.26 


1/18.5 


0.0723346 


0.131552 


0.0065 A 


83.6418 


1/108 


5.91 


1/15.3? 


0.05296 


0.10862 


0.007 A 


90.0757 


1/130 


5.95 


1/18.5 


0.0395829 


0.0911309 


0.0075 A 


96.5097 


1/115 


7.62 


1/16.5 


0.0303889 


0.0780576 


0.007875 A 


?? 


1/115 


7.62 


1/16 


0.02496 


0.07125 


0.008 A 


102.944 


1/220 


2.07 


1/16.5 


0.0233344 


0.0673922 


0.010 A 


128.68 






?? 


0.0101029 


0.0440862 


Slidding BC 




0.0035 ♦ 


45.0379 


1/43.5 


5.2 


1/15.5 


0.913244 


0.920613 


0.0050 4 


64.3398 


1/74 


5.5 


1/18.5 


0.2067 


0.1411 


0.0055 ♦ 


70.7738 


1/80 


3.12 


1/18.5 


0.162815 


0.09994 


0.0075 A 


96.5097 


1/123 


1.07? 


1/17 


0.0303889 


0.0780576 


p = 2 




0.001625* 


41.8209 


1/41 


2.97 


1/17 


0.92876 


0.9317293 


0.00175* 


45.0379 


1/43 


2.77 


1/13 


0.4751098 


0.567067 


0.002* 


51.4719 


1/60 


2.78 


1/17 


0.39088 


0.30088 


0.0025* 


64.3398 


1/75 


2.7 


1/19 


0.139155 


0.22566 


0.0035* 


90.0757 


1/104 


4.72 


1/19 


0.0383266 


0.10147648 


0.0045* 


115.812 


1/190 


2.78 


1/14 


0.0173359 


0.09374663 



TABLE I: Table showing the first two frequency peaks observed in the numerical simulations. All numerical simulations (a) are 
for a = 8 and a mean density p = J \ip\'^dxdy = 1. ■ represents a situation where the density is more or less homogeneous 
in space, that is a liquid. ♦ are numerical simulations with a sliding boundary conditions and * are for p = 2. The two Legget's 
bounds for the superfluid fraction are also computed for each numerical simulation. 

V. CONCLUSIONS AND PERSPECIVES 

We have presented in this review various aspects of the theory of supersolids. A discussion based on fundamental 
physical concepts has shown that a regular lattice with elastic properties may coexist with a superfluid component. We 
derived from the equation of motion of the full system (elastic lattice + superfluid component) the observed property 
that rotation does induce a superflow not implying the whole mass but that a pressure difference does induce a purely 
elastic response without superflow. In this respect the often neglected question of the boundary conditions play a 
central role in the understanding of the properties of the full system. Later on we discussed mathematical models 
of supersolidity, with various degree of realism as far as the detailed microscopic interactions are taken into account. 
Because of the major difSculty of solving honestly and accurately the Schrodinger equation for a sizable number of 
interacting atoms, even with idealized two-body potentials, it is likely that no other fully microscopic theory will 
be available soon for this class of many body systems. We also addressed the connection between the "ideal model" 
leading to the macroscopic equations and the wide variations in superfluid density from one experiment to the other in 
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nominally similar thcrmodynamical conditions. Almost ten years after the experiments by Chan and collaborators and 
many more years after the theoretical predictions by Onsager and by Leggett, one can say that a fair understanding 
of the physics of supersolidity has been reached and likely will be confirmed by experiments now in design or to come. 

The authors acknowledge Peter Mason and Nestor Sepulveda for fruitful discussions. CJ and SR thank the finan- 
cial support of the Agence Nationale de la Recherche through the grant SYSCOM COSTUME ANR-08-SYSC-004 
(Prance). SR thanks the Fondecyt Grant 1100289 (Chile). 
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This is true at low speed. At larger speeds, the velocity difference between normal and superfluid becomes an additional 
thermodynamic parameter. 

As well-known, it is practically impossible to represent accurately in a computer the wave function of N interacting 
particles as soon as A'^ is bigger than a rather small number. Assume that every particle needs, say 3Q points in its 
position space {Q per coordinate), N particles require Q'^^ points of collocation in the position space. Taking Q = 2, 
which is not very large (this amounts to define a function by its value at two points only on the real line!), one gets a 
enormous number of collocation points as soon as A?^ is larger than a few units: for instance N = 8 and Q = 2 yields 
2^* = 1.68 X 10^, meaning that the value of the wavefunction should be known at this number of points, even though 
the the wave function is still quite poorly known (two points per coordinate) and the number of particles is not terribly 
large (eight). This is the well-known problem of computing by brute force properties of quantum systems with more than 
a few interacting particles. In this respect there is a significant difference between the computation of properties like the 
ground-state energy, never very wrong because of the Rayleigh minimum principle behind it, although it is much harder 
to get good representation of functions, like pair correlation for instance, which carry far more information than a number, 
the ground-state energy, without being related to a variational principle. 
[92] Note that in our "perfect crystal" there is one atom per lattice site in the ground state, but that is true only in the limit 

A infinite. Otherwise the assumption "one atom per lattice site" is ambiguous because of quantum effects, see below. 
[93] The potential W{r) is considered as fixed in a first approximation, the atoms others than 1 and 2 staying at rest whilst 
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the atoms 1 and 2 are doing their tunneling trajectories. This is a first order approach to the tunnehng problem, since 
one expects that the other particles do move less and less during the tunneling of 1 and 2 when they are farther and 
farther away. Including this motion of the other particles does not bring any fundamental difficulty, it just makes things 
more cumbersome. The other particles should follow each a homoclinic orbit by starting from and coming back to the 
same site whilst the pair (1, 2) does an exchange trajectory. The tunneling region is everywhere outside of the bottom of 
the potential wells. 

This is shown as follows: if there is no interaction between the two particles they run independent straight trajectories 
moving from one equilibrium to the other and crossing in the middle. Once the interaction V{ri — r2) is turned on, one 
considers all trajectories leaving simultaneously the two equilibria. They are indexed by an angle, if this angle is zero or 
close to zero, the two particles will meet in their mid-course {i.e. near the saddle of the external potential) and fall on 
each other because of the strong attraction at short distance. If on the contrary the angle is large, each particle will fall 
down to infinite depth of the potential well {—AW) without interacting with the other one. In-between, there should be 
a value of the angle such that the two particles make it to the other equilibrium point. 

The function tp{ri,r2) yields in a consistent way the full wave-function of the ground state. One might wonder if, in the 
full crystal, one should not take also into account the tunneling effects between three particles in neighboring sites, etc. 
Although exchange of many particles yields certainly larger actions (and small contributions to the density) exchanges 
between few particles (larger than 2) may contribute significantly to g"" . We plan to investigate this. 
Note that with the standard Lennard-Jones parameters of Helium , ct and e, one has A = (cr/a)^^ ~ 0.6. The two 

body V{r) is needed at shorter distances than usual in physical applications, a fraction of the minimum of a Lennard-Jones 
potential. The l/r^^ repulsion is too strong there because at close distances the dominant interaction is the Coulombian 
repulsion between the nuclei, far less singular than This could affect significantly the estimates of the tunneling 

contribution to the superfiuid density. 

This kind of argument can be extended to show a fluctuation theorem relating the coefficients of response of the super- 
solid, including the superfluid density matrix gf^ to fluctuations in the ground state, via the formal quantization of the 
macroscopic equations of motion. This can be done by using the fact that the phase is the conjugate of the number of 
particles, giving it a well defined meaning. To define also precisely what is meant by the lattice fluctuations, one has deflne 
what is meant by lattice ordering in the ground-state, which can be done by using a method similar to the one leading to 
the supersolid equation in the mean field limit, as done in section [rV] The starting point of the derivation of the quantum 

2 2 

fluctuation theorem is in simple algebraic relations for the quantum harmonic oscillator. Let H — + be the energy 
operator of this oscillator, with a and h positive parameters and p and q non-commuting operators such that [p, q\ — ih. 
Let furthermore (G)q be the mean value of any observable G{p, q) in the ground-state. The eigenvalue of H in the ground- 
state is eo = well-known there is equipartition of energy in this state, so that ^'^V^ = {''^) ~ ^ 

relation allowing to compute a and b as functions of (p^)q and {q'^)^, which yields a kind of fluctuation theorem for the 
harmonic oscillator. To apply this to flnd the response functions of the supersolid (i.e. the Lame coefficients and Qi^), one 
decomposes the possible fluctuations of the macroscopic equations for the supersolid (see section IV I in normal modes. 
This yields a set of harmonic oscillators, four for each wavenumber (four = three elastic modes -I- one phase mode, all 
being linearly coupled), a set to which one can apply the fluctuation "theorem" just derived. The final result is rather 
complex, because of the coupling effects between the different modes. 

Recently P.W. Anderson [Science 324, 631 (2009)] suggested a Gross-Pitaevskii treatment of supersolidity, this description 
is in a completely different frame. Anderson argues, after Andreev and Lifshitz, that, in solid Helium, vacancies undergo a 
Bose-Einstein transition and maybe considered as a dilute (because of the smallness of the superfluid fraction) interacting 
gas which is known to be ruled by the Gross-Pitaevskii equation. Anderson assumes a repulsive interaction, which is true a 
short distance, and neglects the attractive long range elastic interactions between vacancies without rational justiflcation. 
This is in sharp disagreement with the results of Montecarlo simulations showing that vacancies tend to aggregate in 
bubbles [76]. In this respect the "ground-state" proposed by Anderson, assumed to be homogeneous in space, would be 
unstable consisting to an unstable condensate of vacancies. Indeed, Anderson dealing of the condensate of vacancies leads 
to the /ocjismg-nonlinear Schodinger equation. Is well known that this ill-posed equation develops massless flnite time 
singularities. Physically it means that the ground state is made of two phases: a zero density phase with droplets of a 
high density phase where the growth of the long range instability is stopped by the short range repulsive forces, a state 
that is not describable by a weak interaction or mean fleld theory like the non linear Schrodinger equation. 
Note that a Galilean boost with a speed v changes the energy H' = H + P-v+^ mNv^ and the momentum P' — P + mNv 
as usual in classical mechanics. 

Similarly, we may say that by increasing A there is a critical value, Aci , etc. 

All the estimations concern Helium, the energies being in Kelvin {K). It is useful to note that ^ = 12 , the quantum 
of circulation is ^ = 158 m/sA and a — 2.57 A. 

A different trial function that does not satisfy the equation —^f}"{x) — (^)^'i/' could be used: ip{x) = ^1 — (|) 

in x £ [—5, 5] and zero elsewhere. It provides however similar asymptotic results. 
This idea was suggested by Bernard Castaing to Sergio Rica in 1993. 



